OOC  FILE^TO  ADA053169 


THE  COMPUTATION  OF  ECONOMIC  EQUILIBRIA  BY  PATH  METHODS 


BY 


THOMAS  R.  ELKEN 


THE  COMPUTATION  OF  ECONOMIC  EQUILIBRIA  BY  PATH  METHODS  * 


SYSTEMS  OPTIMIZATION  LABORATORY 
DEPARTMENT  OF  OPERATIONS  RESEARCH 


Stanford  University 
Stanford,  California 


Research  and  reproduction  of  this  report  were  partially  supported  by  the  U.S. 

Development  Administration  Contract  EY-76-0B-ni?A  PA 
the  Office  of  Naval  Research  Contracts /NOOOl A- 75-C-G26Z/ aiM  N00014-75-C-j0865l 
the  National  Science  Foundation  Grants*MCS )6-81255  an’u  MCS /6-20'(5T9'. 


Reproduction  in  whole  or  in  part  is  permitted  for  any  purposes  of  the  United 
States  Government.  This  document  has  been  approved  for  public  release  and 
sale;  its  distribution  is  unlimited. 


76^ 


i 


(‘I  // 


I 


CHAPTER 


TABLE  OF  CONTEOTS 


NOTATION  AND  NUMBERING 


INTRODUCTION  AND  SUMMARY 


COMPUTING  ECONOMIC  EQUILIBRIA 


11. 1.  Introduction  

11. 2.  The  Piecewise  Linear  Model  of  Exchange. 

11. 3.  The  General  Piecewise  Linear  Economy... 


THE  BILINEAR  COMPLEMENTARY  PROBLEM  AND  AN 
ALGORITHM  FOR  COMPUTING  EQUILIBRIA  


III.l. 
III. 3. 

111. 3. 

111. 4. 

111. 5. 


The  Bilinear  Complementary  Problem 

The  Algorithm  

Exploiting  the  Linear  Structure 
The  Path -Following  Subroutine  .... 
Moving  From  Cell  to  Cell  


A HOMOTOFY  METHOD  FOR  COMPUTING  BQUILIBKiA. 


IV. 1. 

rv.2. 

IV. 5. 
IV.  4. 


A Convergence  Theorem  

An  Algorithm  Implementing  the  Horaotopy 

Retraction  Method  

The  Path-Following  Algorithm  

The  Cell  Switching  Decisions  and 
Operations  


RESULTS  FROM  COMPUTATIONAL  EXPERIMENTS 


Computer  Implementation  of  the  Algorithms. 

Experimental  Design  

Numerical  Results  

Conclusions  


SUGGESTIONS  FOR  FUTURE  RESEARCH 


APPENDIX:  CONSTRUCTING  A PIECEWISE  LINEAR 
APPROXIMATION  


BIBLIOGRAPHY 


ii 


NOTATION  AND  NU^^BERING 


All 

1) 

2) 

3) 

h) 


5) 


6) 

7) 


8) 


9) 


We  define  some  notation  which  may  not  be  completely  standard, 
points  and  sets  are  in  unless  indicated  otherwise. 

b"  H {x  e e”|x  > 0] 

n 

For  X,  y £ 3R  , denote  their  inner  product  by  (x,y)  = Z x.y. . 

i=l  ^ ^ 

1 /p 

Let  I|x||  = (x,x)  ' be  the  norm  of  x. 

Denote  by  d(*,-)  the  distance  function  defined  as 

d(A,B)  = inf  ||x-yl| 
x£A 
y€B 

Let  the  ball  about  A of  radius  e , where  A may  be  a point  or 
a set,  be 

B(A,e)  = {x  e IR’^|d(x,A)  < e] 

For  any  positive  integer  m,  let  m=  {1,2,  ...,m).  If  m = 0, 
then  m s 

conv[A^,A^, . . . ,A^]  is  the  convex  hull  of  the  sets  (or  points) 

A^,  A^,  ...  , A^. 

j^nxm  space  of  all  real-valued  matrices  with  n rows 

and  m columns. 

If  a c n and  p cz  m are  index  sets,  then,  for  A G IR 


iii 


consists  of  the  rows  of  A indexed  by  a,  and 


Numbering 

The  chapters  are  numbered  by  Roman  numerals  and  the  sections 
are  numbered  consecutively  within  each  chapter,  i.e.,  II. 1,  11,2,  etc. 
All  theorems,  lemmas,  and  examples  are  numbered  consecutively  within 
each  section.  Equations  are  numbered  separately  and  are  identified 
by  being  enclosed  in  parentheses.  Equation  (V.1.2)  is  the  second 
equation  in  Section  V.l,  for  instance.  The  chapter  numeral  is  omitted 
for  results  or  equations  referred  to  within  the  same  chapter. 
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Thomas  R.  Elken 

ABSTRACT 

An  introduction  to  the  economic  equilibrium  model  is  given 
and  it  is  demonstrated  that  a path  method  can  be  used  to  compute 
equilibria  for  pure  exchange  economies  in  a nonlinear  setting. 

Next,  a model  is  described  for  an  economy  in  which  the  utility 
functions  are  piecewise  linear  and  the  cons’jmption  and  production  sets 
are  polyhedral.  It  is  shown  that  an  equilibrium  for  this  economy  is 
the  solution  to  a system  of  bilinear  equations  subject  to  certain 
linear  inequality  and  complementarity  constraints. 

Two  approaches  are  discussed  for  computing  equilibria  for  such 
economies.  The  first  is  the  bilinear  complementarity  algorithm  (BCA) 
and  the  second  is  the  homotopy  retraction  algorithm  (HRA).  Convergence 
proofs  are  given  for  both  methods  using  the  general  theory  for  path 
methods  described  above. 

The  BCA  and  HRA  have  been  implemented  as  computer  progratr^. 
Detailed  descriptions  of  the  algorithms  are  given,  and  the  results  of 
some  numerical  experiments  are  reported.  Seven  small  problems  were 
solved  by  both  algorithms.  No  conclusion  could  be  drawn  as  to  which 
algorithm  was  superior,  but  both  performed  well  enough  that  it  appears 
that  much  larger  equilibrium  problems  also  can  be  solved  efficiently 


by  these  methods. 


THE  COMPUTATION  OF  ECONOMIC  EQUILIBRIA  BY  PATH  METHODS 


CHAPTER  I 

I.  Introduction  and  Summary 

This  report  is  concerned  with  two  algorithms  for  computing 
equilibria  for  a general  piecewise  linear  economy.  Proofs  will  be  pre- 
sented which  show  that  these  algorithms  converge  under  certain  conditions 
and  the  computational  results  from  the  first  implementation  of  these 
algorithms  will  be  reported.  The  first  algorithm  is  the  bilinear  com- 
plementarity algorithm  of  Wilson  [1976].  A new  proof  of  convergence 
is  given  here  using  the  results  on  subdivided  complexes  in  Elken  [1977], 
and  a detailed  description  of  an  algorithm  which  is  suitable  for  numerical 
implementation  is  presented.  The  second  algorithm  is  new  but  it  is 
heavily  influenced  by  the  ideas  of  Wilson  [1976]  and  Kellogg,  Li,  and 
Yorke  [1976].  We  call  it  the  homotopy  retraction  algorithm.  Both 
algorithms  show  great  promise  in  computing  equilibria  for  larger  models 
than  has  been  possible  up  to  now. 

This  report  constitutes  the  second  half  of  the  author' s dissertation. 
The  first  half  is  contained  in  Elken  [1977];  we  will  call  this  work 
Part  1.  The  reader  who  is  interested  in  the  proofs  for  the  two  convergence 
theorems  in  this  report  will  have  to  be  familiar  with  the  definitions 
and  results  in  Chapter  II  of  Elken  [1977].  Those  readers  interested  in 
the  algorithms,  their  implementation,  and  the  numerical  experiments  will 
find  that  this  report  contains  all  the  relevant  material. 
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In  Chapter  II  of  this  report  we  Introduce  the  economic  equilibrium 


problem  in  its  traditional  form  and  show  how  it  can  be  solved  by  a path 
method.  (A  path  method  is  a procedure  for  solving  a system  of  nonlinear 
equations  by  following  a (piecewise)  differentiable  path  from  a known 
starting  point  to  a solution  for  the  system  of  equations;  see  Elken  [1977] 
for  more  details.)  We  also  introduce  the  piecewise  linear  formulation 
of  an  economy  as  it  was  developed  by  R.  Mantel  [I967]  and  R.  Wilson  [1976]. 
It  is  this  formulation  which  we  use  to  develop  algorithms  to  exploit  the 
linear  structui'e  to  compute  equilibria  for  this  model  of  an  economy. 

In  Chapter  III  we  present  a generalization  of  the  piecewise  linear 
equilibrium  problem,  the  bilinesu:  complementarity  problem  (BCP)  of  Wilson 
[1976].  Wilson  defined  a path  and  proved  that  it  led  from  an  easily 
obtained  starting  point  to  a solution  of  the  equilibrium  problem.  We 
present  a new  proof  of  this  and  an  algorithm  which  exploits  the  linearity 
of  this  model  to  reduce  the  nonlinear  path-following  problem  to  the 
lowest  possible  dimension  (<  the  number  of  consumers) . The  last  part  of 
this  chapter  is  a description  of  the  algorithm  as  it  is  implemented 
in  a computer  program. 

In  Chapter  IV  we  present  a new  path  method  for  computing  economic 
equilibria  and  prove  that  it  is  convergent.  We  also  describe,  in  detail, 
an  algorithm  based  upon  this  path  method  which  has  been  numerically 
tested. 

Chapter  V contains  the  results  of  some  computational  experiments 
with  the  algorithms  presented  in  the  previous  two  chapters.  The 
conclusions  are  primarily  a comparison  of  these  two  methods,  but  it  is 
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hoped  that  comparisons  with  existing  codes  be  made  in  the 

future. 

Chapter  VI  suggests  some  promising  directions  for  future  research. 
An  Appendix  is  included  which  describes  one  method  for  generating  piece- 
wise  linear  approximations  to  nonlinear  utilit;y  functions  of  several 
variables. 
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CHAPTER  II 


t 


COMPUTING  ECONOMIC  EQUILIBRIA 


II. 1.  Introduction 

The  concept  "equilibrium"  suggests  that  opposing  forces  are 
in  balance.  Economic  equilibrium  theory  is  concerned  with  the  balance 
of  the  demand  from  consumers  with  the  supply  provided  by  producers. 

Each  agent  in  the  economy  is  concerned  only  with  the  maximization  of 
his  satisfaction  or  profit  depending  upon  whether  he  is  a consumer  or 
producer,  respectively.  To  be  more  explicit,  a special  class  of  economies 
is  considered,  namely,  private  ownership  economies  in  which  consumers  own 
the  resources  and  control  the  producers.  Given  a price  system,  each 
producer  maximizes  his  profit,  which  is  distributed  to  cons\imer- share- 
holders. The  wealths  of  the  latter  are  thereby  determined,  and  they 
maximize  their  satisfaction  frcaa  the  consun5)tion  of  goods  subject  to 
their  wealth  constraints.  As  a result  of  this  process,  each  agent  chooses 
an  action.  The  state  of  an  economy  in  which  these  optimal  actions  of 
consumers  and  producers  are  compatible  with  the  resources  is  called 
an  equilibrium. 

The  traditional  questions  ^ich  have  been  studied  by  economists 
are:  Does  a price  system  exist  which  puts  the  economy  in  equilibrium? 

Is  an  economy  at  equilibrium  stable?  That  is,  if  there  is  a violation 
of  one  of  the  conditions  of  equilibrium,  do  rational  actions  of  the 
agents  tend  to  restore  an  equilibriian.  An  excellent  work  which  deals 
with  these  questions  is  Arrow  and  Hahn's  General  Congetltive  Analysis 
[ 197l] . This  work  also  provides  an  excellent  introduction  to  the 
literature  of  economic  equilibria  theory. 
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This  work  is  concerned  with  the  existence  question  as  it 


arises  in  the  prohlera  of  computing  equilibri’jm  prices  and  activities. 

The  question  of  stability  will,  not  be  discussed.  The  ability  to  compute 

equilibria  for  large-scale  models  is  important  because  it  may  allow  a 

great  improvement  in  economic  modeling. 

The  difficulties  inherent  in  economic  model  are  well  known. 

"The  problem  of  collecting  reliable  data  on  the  technological 
processes  that  are  currently  available  is  enormous,  to  say 
nothing  of  the  difficulty  of  inventing  appropriate  input- 
output  coefficients  for  productive  techniques  that  remain 
to  be  discovered.”  (Scarf  [1973],  pp.  7-8). 

Because  of  the  large  size  of  most  economic  models,  the  linear  programming 
formulation  is  invariably  used.  The  advantage  of  the  equilibrium  model 
over  linear  programming  is  that  consumer  demands  and  their  dependence 
upon  price  are  recognized  and  modeled  with  some  rationality.  Unless 
the  system  is  prepared  to  tolera-te  strict  rationing,  considers  will 
respond  not  only  to  the  availability  of  items  but  also  to  their  price. 

Thus,  we  see  the  importance  of  computing  equilibria  as  an  aid  in 
economic  planning. 

The  first  algorithms  for  the  computation  of  economic  equilibria 
were  developed  by  Scarf  [ 1967 ] • He  developed  an  algorithm  for  solving 
fixed  points  of  a function  which  maps  the  unit  simplex  into  itself. 

Then  the  equilibrium  problem  was  transformed  into  such  a fixed  point 

j 

problem.  To  illustrate  how  this  transformation  can  be  accomplished  j 

we  will  deal  with  a pure  exchange  economy,  i.e.,  we  ignore  production  | 

for  the  moment.  I 
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Suppose  there  are  n goods  in  the  economy  of  m consumers 

or  traders,  and  that  each  of  the  consumer's  preferences  are  represented 

by  a utility  function.  To  be  precise,  let  y e be  a vector  of  the 

n goods  and  tt  £ IR  be  the  vector  of  prices  for  the  n goods.  A 

bundle  of  goods  x is  preferred  to  a bundle  y by  consumer  i, 

i = 1,  ...,m  if  and  only  if  u^(x)  > u^(y)  where  ® is  a 

strictly  concave  and  continuous  utility  function.  The  demands  of  the 
"til 

i consumer  are  determined  by  the  solution  to  the  following  problem; 

maximize  u^(y) 

subject  to  Try  < ttw^  , (l.O) 

where  w is  the  initial  endowment  of  the  i consumer,  i = 1,  ..,,m. 

We  shall  assume  that  the  solution  to  this  problem  can  be  written  as  a 
continuous  function  of  the  prices  tt,  d^(7r).  The  individual  trader's 
excess  demand  function  is  d^(Tr)  - w^,  i = m.  The  excess  demand 

will  be  positive  for  those  commodities  whose  stock  he  wishes  to  increase 
by  exchange  and  negative  for  the  remaining  items.  The  market  excess 
demand  function  g is  the  sum  of  the  individual  excess  demand  functions 

gCir)  = E (d  (tt)  - w^) 
i=l 

An  equilibrium  price  vector  tt  is  one  for  which  all  of  the 
market  excess  demands  are  less  than  or  equal  to  zero 

g(7r)  < 0,  (l.l) 

TT  > 0,  TT  5^  0 , (1.2) 

T 

TT  gCir)  = 0 . 
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(1.5) 


Equation  (1.3)  implies  that  there  is  a zero  price  for  any  commodity 
whose  demand  is  strictly  less  than  zero. 

An  important  property  of  the  excess  demand  function  is  that  (l.3) 
will  be  satisfied  for  any  price  system  v because  the  value  of  each 
consximer's  excess  demand  is  identically  zero.  This  property  is  known 
as  Walras'  law.  Another  property  which  is  easy  to  derive  is  that  the 
excess  demand  function  is  homogeneous  of  degree  zero,  that  is,  if  the 
prices  are  scaled  up  by  a constant  factor,  the  demands  remain  the  same. 
This  implies  that  it  is  sufficient  to  search  for  an  equilibrium  price 
vector  on  the  unit  simplex,  S^.  The  following  function  from  into 

has  the  property  that  a fixed  point  is  an  equilibrium; 

X + max(0,g  (tt)) 

= , i t n . 

1 + E max(0,g  (tt)  ) 
i=l 

To  see  that  a fixed  point  of  f is  an  equilibrium  point, 
the  equation  tt  = f(7r)  can  be  written 

OTr.  = TT.  + max[0,  g.  (tt)  1 

with  a = 1 + E^  max[0,g.  (tt)  ].  If  a is  in  fact  greater  than  one, 
then  the  condition  Tr^(a-l)  = max[0,g^(Tr)  ] implies  that  g^(7r)  > 0 
vdienever  > 0.  Since  some  tt^  are  strictly  positive  this  violates 
TT  g(v)  = 0.  Therefore  a = 1 and,  hence,  g^(u)  < 0 for  all  i, 
and  TT  is  an  equilibrium  vector. 


One  procedure  for  moving  towards  equilibrium  is  a price  adjust- 


ment in  which  the  price  of  a good  is  increased  if  the  excess  demand 
for  that  good  is  positive,  decreased  if  the  excess  demand  is  negative. 

This  is  the  classical  tatonnement  or  "groping"  for  an  equilibrium 
(Walras  [1874]).  Scarf  [i960]  has  shown  that  for  virtually  arbitrary 
excess  demand  functions,  this  process  can  be  unstable.  This  price 
adjustment  is  globally  stable,  however,  if  a certain  gross  substitu- 
tability between  all  commodities  is  satisfied  (Arrow  and  Hurwicz  [1958]). 

The  differential  equation  which  expresses  this  process  is 

||  = sM  (1.4) 

We  vri.ll  ncv.'  show  that  a convergent  price  adjustment  process  can 
be  defined  which  behaves  like  (1.4)  initially.  Suppose  that  the  excess 
demand  function  g is  differentiable.  Since  the  excess  demand  function 
is  homogeneous  of  degree  zero,  g' (tt)  has  rank  n-1,  in  general;  hence, 
we  must  re-^uce  the  dimension  of  the  problem  by  one.  This  can  be  done 
by  solving  the  problem  on  the  unit  simplex  or  the  nonnegative  portion 
of  the  unit  sphere,  [x  € ]R^|x  > 0,  ||x||  =1).  We  shall  follow  a certain 
amount  of  tradition  in  mathematical  economics  by  considering  a distinguished 
commodity  or  'n’jmeraire'  (Quirk  and  Saposnik  [I968])  so  that  the  price 
of  other  goods  are  measured  in  terms  of  this  numeraire  good. 

Suppose  there  are  n+1  commodities  -with  prices 
vdiere  tTq  is  the  unit  price  of  the  numeraire  good.  We  make  the  follow- 
ing assumption  on  the  excess  demand  function  g^(7r),  i = 0,  ...,n. 
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Asstimption  1.1.  If  tt^  =0,  then  g^iv)  >0,  i = 0, 1, 

Thus,  no  equilibrium  price  vector  will  be  on  the  boundary  of 
So  we  will  work  in  the  domain  of  price  systems  {tTq,  , . . ,tt^) 
which  satisfy  tTq  > 0.  In  this  domain,  by  the  homogeneity  property 
of  price  systems,  one  can  normalize  a price  vector  by  dividing  by  tTq. 

Thus  a price  system  (ttq,  . • . ,Tr^)  can  be  represented  by  a xinique  point 

(Pl>---/Pji)  = •••  > ® + 

With  this  interpretation,  the  space  of  price  systems  is 

2 

Suppose  that  the  excess  demand  functions  g^Cir),  ...  , g^Cir)  are  C 
and  satisfy  Walras'  law.  Then 

gQ(7r)  = - E g.  (tt)  . (1.5) 

i=l  ^0 

This  along  with  the  boixndary  condition  1.1  implies  that  the  equations 
g^M  = 0 , i = 1, n 

are  necessary  and  sufficient  conditions  for  economic  equilibrium. 

Let  fj^(p)  = g^(l,  tt^/tTq,  ...  , tTjj/tTq)  , i = 1,  ...,n  . 

Then  consider  the  application  of  the  following  version  of  Kellogg, 
Li,  and  Yorke's  deformation.  Pick  some  p*^  such  that  p^  € ]R^  and 
for  exactly  one  i £ n,  p^  = 0,  say  i = 1.  Then  the  boxmdary  condition 
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(III. 1.1)  can  be  shown  to  hold  on  D = {p|o  < < P,  i 6 n)  for  some 

P > 0 because  of  Assumption  1.1  and  equation  (1.5)  (cf.  Smale  [1976], 
p.  Il6).  M = D X [-1, 1]  and  is  denoted  by 

F(p,0)  = 9f(p)  - (l-e)(p-P°)  . 


Then  by  the  theory  of  Section  II. 1,  if  0 is  a good  value  for  F, 
there  is  a curve  with  boundary  points  at  (p^,0)  and  (p*, 1) 

and,  of  course,  f(p*)  = 0.  But  if  (p(t),e(t))  = [0,T]  -» F ^(O)  such 
that  (p(0),e(0))  = (p°,0)  and  (p(T),0(T))  = (p*,l),  then  i)(t) 


satisfies 


or 


= 0 


(tf’  (p)  - (l-0)llf(p)  - (p-p  )] 


p(t)' 

e(t) 


= 0 . 


At  t = 0 we  have  0=0,  p = p , 5(t)  = A > 0,  and 


p(t)  = (p) 


(1.6) 


Thus,  initially  at  least,  the  adjustment  of  prices  by  a retraction 
method,  is  the  same  as  the  classical  tatonnement  method  of  price 
adjustments.  It  can  be  said,  with  reference  to  the  equilibrium  problem, 
that  the  parameter  0 of  the  retraction  method  allows  a smooth  transition 
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[0f’(p)  - (i-y)i]p(t)  = -a(t)(f(p)  - p + p°)  , 
or 

f (p)  p(t)  - p(t)  = - (f(p)  - p + p°)  (1.7) 

As  p(t)  -> p*  and  0(t)  -> 1 along  F ^(O)  it  is  clear  that  the 
second  term  on  the  left  side  of  (5.7)  vanishes.  Also  for  (p, 0)  £ F ^(O), 
p - p®  = 0/(i-0)  f(p),  hence,  (5.7)  is  approximately 

f’(p)  P(t)  = - 0(t)  (|  + ~)  f(p)  . 

So  for  some  'K>0,  p(t)  almost  satisfies  f (p)  f)(t)  = -Xf(p),  the 
global  Newton  equation  of  Smale  [ 1976] . 

II. 2.  The  Piecewise-Llnear  Model  of  Exchange 

The  usual  assun^itions  made  in  the  literature  of  economic  equi- 
librium theory  include  asstia^tions  that  the  utility  functions  are 
concave,  the  consumption  sets  are  convex  as  are  the  production  possi- 
bility sets.  It  is  well-known  that  concave  functions  can  be  approximated 
to  an  arbitrarily  small  tolerance  by  piecewise  linear  functions  and 
that  convex  sets  can  be  approximated  very  closely  by  polyhedral  convex 


from  tatonnement  to  a global  Newton  method.  To  see  this  consider  eigain 
the  relation 


! 
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sets.  In  this  section  the  discussion  will  center  around  a model  of  an 
economy  in  which  the  utility  functions  are  assumed  to  be  concave  and 
piecewise  linear,  and  the  consumption  and  production  sets  are  polyhedral. 

This  model  was  introduced  by  Rolf  Mantel  [ I968] . He  gave  an 
ingenious  proof  of  the  existence  of  economic  equilibria  without  using 
the  methods  of  combinatorial  or  differential  topology.  Mantel' s 
approach  did  not  seem  to  be  computationally  efficient,  for  it  involved 
a con^Jlex  limiting  operation  at  one  point.  G.  Dantzig,  B.  C.  Eaves, 
auid  D.  Gale  [1976]  did  use  this  model  as  the  basis  for  a new  approach 
to  computing  equilibria.  They  solve  the  problem  by  computing  a fixed 
point  of  a point-to-set  map  whose  values  are  determined  by  the  solution 
to  a linear  program.  We  will  not  discuss  this  algorithm  further, 
even  though  it  promises  to  be  one  of  the  main  competitors  to  the 
algorithms  we  present  here. 

Both  Dantzig,  Eaves  and  Gale  [1976]  and  this  work  are  concerned 
with  solving  equilibrium  problems  with  a relatively  small  number  of 
households  or  traders  (about  3-10)  and  a large  number  of  goods  (up  to 
300).  One  obvious  example  of  a problem  of  such  a scale  would  be  if  3 
to  10  co\mtries  were  considered  to  be  consumers  who  were  involved  in 
the  production  and  trade  of  up  to  300  goods. 

For  simplicity,  we  will  first  present  a model  which  allows  only 
linear  utility  functions  and  ignores  production.  Later  it  will  be 
shown  how  the  general  piecewise-linear  economy  can  be  formulated  and 
solved  in  essentially  the  same  manner. 
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Consider  a simple  exchange  economy  with  m households  and 


£ goods.  Each  household  i £ m has  available  a finite  set  S.  of 

i ^i 

activities  such  that,  if  it  chooses  a vector  z £ K of  nonnegative 
activity  levels,  then  it  obtains  the  utility  and  it  consumes 

the  vector  B z of  quantities  of  the  commodities,  where  y £ K ^ 
and  If  w^  £ E is  household,  i's  initial  endowment 

of  goods,  then  its  excess  demand  function  is  g^(z^)  = B^z^  - w^. 

Notice  that  in  the  terminology  of  the  pure  exchange  problem  in  Section  1, 
U^(z^)  = y^z^.  Also,  in  that  pxire  exchange  problem,  if  all  of  the 
activities  in  are  merely  the  consumption  of  a particular  commodity 

then  all  of  the  technology  matrices  are  identities  (B  = I £ E ) . 

Definition  2.1.  A price  vector  tt  ^ 0,  and  a consumption  allocation, 
z^,  i £ m constitute  a competitive  equilibrium  iff 

a)  the  net  trades  are  feasible  with  free  disposal: 

ZBz<Sw.,  z‘’'>0,  i£n^ 
i=l  i=l  ^ 

and 

b)  Given  the  nonnegative  price  vector  tt,  z^ 

maximizes  U(z^)  = y^z^ 
subject  to  ttB^z^  < ttw^  , z^  > 0 , 

for  each  i £ m. 

This  notion  of  equilibrium  is  clearly  equivalent  to  that  in 
Section  1.  If  z^  is  considered  to  be  a function  of  tt,  then  B^z^ 


A similar  notion  is  that  of  a quasi-equilihrium,  vrtiich  is 


defined  as 

Definition  2.2.  A specification  (it,  z^,  i G m)  is  a quasi-equilibrium 
iff 

m . . m . 

a)  Z < E w^,  z > 0,  i e m and 

i=l  i=l 

b)  Given 

minimizes  7r(B^z^  - w^) 

^ i a.  i i ^ i ^ ^ 

subject  to  yz  >yz,z  >0 

for  each  i € m,  and  this  minimum  is  zero. 

Condition  b)  can  be  interpreted  as  choosing  z^  to  minimize 
the  net  expense  of  maintaining  the  utility  level  U^(z^)  and  requiring 
that  the  budget  be  balanced  exactly  7r(B^z^  - w^)  =0. 

The  following  conditions  will  be  assumed  throughout* 

Assumption  2.3. 

a)  For  each  consumer  i £ ^ the  consvutption  set 

= {x  £ IR"^!  (3  z^  > 0)  B^z^  < x] 

is  bounded  below. 

b)  The  induced  utility  function  U^(x)  = maxly^z^j  z^  >0,  B^z^<x]  on 

is  insatiable,  i.e.,  for  all  x^  £ X^^  3 x^  £ X^  such  that 
U^(x^)  < U^(x^),  i G m. 

c)  The  initial  endowment  w^  is  strictly  positive  for  each  i £ ra. 

1I4 

J 


Debreu  [1959]  has  shown  that  these  assumptions  imply  that  a 
competitive  eq\iilihrium  exists,  and  that  a quasi-equilibrium  is  equivalent 


to  an  equilibria. 

Let  us  consider,  for  a moment  the  complexity  of  the  problem  of 
computing  equilibria  for  this  class  of  economies.  B.  C.  Eaves  [1975] 
has  shown  that  if  = I,  i £ ^ this  problem  can  be  formulated  as  a 
linear  complementarity  problem,  and  with  these  assumptions.  Lemke' s 
algorithm  [19^5]  will  yield  a solution  after  a finite  n\amber  of 
additions,  multiplications  and  comparisons.  However,  when  B^  is 
allowed  to  be  more  general,  a problem  with  rational  data  may  have  an 
irrational  solution.  The  following  example,  due  to  Andreu  Mas-Colell, 
is  such  a problem: 


= (-25) 


i 

w = 


i = 1,2,3. 


Mas-Colell  has  shown  that  every  equilibrium  price  vector  w is  a 
positive  scale  of  (l  + y/J,  1).  Thus,  no  finite  algorithm  can  hope 
to  solve  this  class  of  equilibrium  problems  exactly. 

Next,  we  introduce  the  auxiliary  linear  program  which  has  the 
property  that,  if  the  correct  constants  are  used  in  the  right-hand  side, 
the  solution  to  the  linear  program  is  a quasi-equilibrium  and,  hence, 
an  equilibrium. 
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The  objective  of  the  linear  program  is  to  maximize  p,  a measure 


of  exports  from  the  system,  subject  to  maintaining  utility  levels  v. 


I 


and  feasibility  of  trades,  i.e., 

maximize  p 


subject  to  i £ m 


Z B^z^  + pe  < Z 
i£m  i=l 


m 


(2.1) 


icm 


z > 0,  i £ m , 


vrtiere  e £ E is  a strictly  positive  vector.  This  formulation  is 
due  to  Wilson  [1976].  It  also  was  influenced  by  Debreu's  [1951]  notion 
of  a coefficient  of  resource  utilization.  If  e = Z™  , w^,  then 
cp  = 1-p  is  the  coefficient  of  resource  utilization.  In  Debreu’ s 
theory,  if  cp  = 1,  then  the  economy  is  on  the  Pareto  optimal  frontier 
(any  increase  in  one  consumer's  utility  would  require  another  consumer's 
utility  to  decrease  for  the  trades  to  remain  feasible) . 

The  d\ial  for  the  auxiliary  linear  program  is 


Here,  v is  the  price  vector  and  can  be  interpreted  as  i's 

marginal  utility  of  income.  Let  u^  s 7)(w^  - i £ be  i's 

budget  stirplus.  If  compieraentajry  slackness  holds  for  the  dual  problems 
(2.1),  (2.2)  then  we  have  u^  = ttw^  - A^v^ 

for  i £ m.  Also, 

7r(E  B^z^  + pe)  = 7r(Z  w^) 
i i 


and 


imply 


ire  = 1 


p = E u^  • 

i 


(2.3) 


Thus,  at  a solution  to  the  primal  and  dual  programs,  these 
relationships  will  hold. 

Lemma  2.4.  If  the  utility  levels  v^,  are  chosen  so  that  u^  = 0,  i £ ^ 
at  an  optimal  solution  (v;  i £ ra;  z^,  i £ m,  p)  to  the  problems 
(2.1),  (2.2),  then  (tt;  z^,  i £ m)  is  a quasi -equilibrium. 


Proof.  Clearly,  the  trades  are  feasible  because  (2.3)  implies  p = 0. 

All  that  remains  to  show  is  that  z^  minimizes  7r(B^z^  - w^) 

subject  to  r^z~  > y^z^  = v^^,  z^  > 0.  Call  this  problem  (P) . The  dual 

i ••  i “•! 

problem  (D)  is  to  maximize  subject  to  , A^  > 0.  z 

and  Aj,  satisfy  the  constraints  of  (P)  and  (D)  because  they  satisfy 
the  constraints  for  (2.1)  and  (2.2).  Also,  the  complementary  slackness 
conditions,  A^^r^z^  = ^v^  and  = irB^z^  are  satisfied  by  the 

optimality  of  the  variables  for  (2.1),  (2.2).  Thus,  z^  is  an  optimal 
solution  of  (P).  I 


17 


This  lemma  provides  motivation  to  parametrically  vary  the 

initial  utility  assignments  v^  and  adjust  the  solutions  of  (2.1) 

and  (2.2)  xintil  u^  = 0 for  all  i € m.  Consider  slack  variables 

t^  > 0 such  that  “ ^i  ^ ^i"  instead  of  varying  v^, 

we  cam  vary  t^,  ignoring  the  complementairy  slackness  with  the  dual 

variable  X. . 

1 

Thus,  we  can  state  the  equilibria  problem  as: 

Find  (i  G p;  s,  t,  (i  € m),p)  which  satisfy 

i i . . ^ 

y z -t.  =v.  itm, 

11  — 

m . . -t 

ESz^  + pe  +s=Ew, 

i=l  i=l 

ttB^  - A^r^-  = 0 , i e ^ 

Tre  - p = 1 , 

TTS  = 0 , 

^^z^  = 0 , i € m , 

and  (tt,  (i  £ m),  pj  s,  t,  z^  (i  £ m,  p)  > 0 


and 


i 

TTW 


Ai(Vi  + t^)  = 0 , 


i £ m , 


(2.4) 


Solution  procedures  for  this  formulation  of  the  problem  will  be 
discussed  below,  but  first,  it  is  important  to  discuss  the  choice  of 
the  initial  utility  levels  v^,  i £ m. 

An  obvious  candidate  for  the  initial  utility  level  assignment 
Vj^  is  consumer  i 's  induced  utility  vf  for  the  initial  endowment  w^, 
i.e.  , 
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vf  = Uj^(w^)  = inax{-ir^z^ I > 0,  B^z^  < w^)  , 


(2.5) 


or  dually 

= min{7rV|Tr^  > 0,  ttV  > (2.6) 

Lemma  2. 5.  With  the  choice  = vf,  > 0,  for  i £ m at  a solution 
to  (2.3). 

Proof,  u^  = TTW^  - by  definition.  Suppose  > ttw^.  Since 

TTw^  > 0,  we  cannot  have  = 0.  If  A^  > 0,  we  have 

V*  > aT^ttw^  . 

1 1 

However,  A^^V  > 0,  and  by  the  dual  feasibility  (2.2)  of  A^  and  ir, 

A^  vB  > r 

which  contradicts  the  optimality  of  v*  in  (2.6).  Hence, 

U.  = TTW^  - A.v?-  > 0.  I 
111—" 

Since  we  are  trying  to  find  a solution  of  (2.5)  for  which  u^  = 0, 
it  would  seem  to  be  advantageous  to  choose  v^  as  large  as  possible. 

For  arbitrary  choices  of  v^^  larger  than  v*  the  constraints  (2.5) 
may  not  be  feasible.  For  theoretical  purposes  it  turns  out  to  be 
preferable  to  choose  v^  G (O^v*).  In  this  case,  the  following 
equivalent  properties  are  easy  to  verify 
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d 


r 


(1)  A.t.  = 0 implies  u.  = ttw^  - A.  (v.  + t . ) >0 

i (2.7) 

(2)  u.  = TTW  - A.  (v.  + t. ) < 0 implies  A.t.  > 0 . 

' ' 1 ill—  11 

This  property  will  be  exploited  from  two  points  of  view; 

(a)  If  the  condition  > 0 is  enforced,  (2.7)  is  a complementary 
type  of  relationship  between  and  i £ 

(2)  If  one  considers  f(A;7r,t)  s t^w^  - + t^^)  and  one  is  solving 

f(Ajr, t)  = 0 subject  to  A £ then  (2.7)  is  the  boundary  con- 
dition that  states,  for  A £ f(X;7r,  t)  points  into  ]R^. 

11,3.  The  General  Piecewise  Linear  Economy 

Now  we  allow  further  constraints  to  be  added  to  each  consumer's 
consumption  set  and  a finite  number  q of  firms  which  are  owned  by 
the  consumers.  Suppose  that  each  consumer  i owns  the  fraction  cr^ 
of  firm  J,  where  cr^  > 0,  J £ q,  and  cr^  = 1,  for  each  J £ q. 

Before  we  can  define  the  consumption  sets  we  must  define  production 
sets  for  each  of  the  q firms.  Let 

y.  = {y  £ nR^|(3  u^'  > 0)  D^u'^  < d^,  y < 

0 

th  i 

be  the  J firms'  production  set  where  or  is  a vector  of  initial 

endowments  for  each  J £ ^.  Each  household's  consumption  set  is  defined 

as 
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X.  = {x£  ]R^|  (3z^  > 0)  < a^,  X > - E aV,  y^'  € Y ) . 

^ “ jea  ^ 


With  these  definitions,  we  can  define  an  equilibrium  for  economies 
with  production 


Definition  3.1.  A price  vector  ir,  a consumption  allocation  z^,  i £ rm 
and  a production  allocation  y^,  J £ £,  constitute  a competitive  equi- 
librium if  and  only  if 


a)  TT  > 0 and  ~r  ^ 0, 


m . . q • m . 

b)  Z < Z y^  + Z W^ 


i=l 


c)  y^  maximizes  Try 

subject  to  y € Y.,  j £ q,  and 
d 

-i  • • i i 

d)  z maximizes  r z 

subject  to  ttB^z^  < 7r(w^  + Z u^y^)  , A^z^  < a^,  for  every  i £ m. 

jsa  ^ 

A quasi-equilibrium  (or  the  compensated  equilibrium  of  Arrow  and  Hahn 
(1971])  replaces  condition  (d)  with 

d' ) z^  minimizes  ttB^z^ 

subject  to  r^z^  > , and  Tr(B^z^  - w^  - Z aty^)  = C,  i £ m. 

jea  ^ 

If  we  combine  Assiimption  2.3  (using  the  new  definition  of  consumption  set) 
with  the  assumption  below,  then  equilibria  exist  and  coincide  with  quasi- 
equilibria (Debreu  [ 1959] ) • 


i=l 


i=l 
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Assumption  3.2. 

a)  free  disposal 

b)  production  is  irreversible 


c)  (Y.  - n]R  ={0),  i.e.,  there  is  no  free  production  and  the 

.1  + 


firm  may  produce  nothing  at  all. 

The  auxiliary  linear  program  is  now  defined  as 


maximize 

P 

subject  to 

i i 

r z 

>Vi  , 

„ i i 

i 

A z 

< a , 

D*^u^ 

< d^. 

m . . <1  4 . m . q . 

Z - S E^u^  + pe  < Z w^  + Z O)^ 

i=l 


i=l  J=1 


> 0,  i £ u^  > 0,  j £ q. 


whose  dual  is 


m 


minimize  L,  [ttw  + ^<^a^  - 6,.v^]  + Z B'^d 
. , h h . , 

1=1  J=1 


subject  to 


me  = 1 


- ttE^  > 0 


+ ttB^  - > 0 


7r,  B^  (j  £ £),  V , A.  (i  £ m)  all  > 0 


Next  we  prove  the 
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(5.1) 


(3.2) 


r 

1 

i 


i 


Lemma  3.3.  At  a solution  to  (5.1)  and  (3.2)  and  for  a solution 
to  max{7ry|y  G Y^},  the  following  equations  holdi 

u.  = 7f[w^  + E cr^y*^  - B^z^] 

= Tr(v^  + t auj^)  + + E 6^(o-^d^)  - A (v.  + t ),  for  all  i € m. 

0=1  ^ 0=1  ^ 111 

Proof,  y^  satisfies  iry^  = max  Try 

subo'ect  to  D^u^'  < d*^  (P) 

y - E"u^  < (xP , u.  > 0 . 

- 0 - 

Hence,  iry^  = min  5^d*^  + ttlo^ 

s.t.S^D^  - ttE^  > 0 (D) 

IT  = rr,  6^'  > 0 . 


But,  using  the  complementary  slackness  properties  of  the  optimal  solution 
(z^,  u^*,  p;  TT,  5^',  v*^)  to  (3.1)  and  (3.2),  it  is  true  that 


6^L^u^  - fE^u^  = 0 


and 


S^D^u^  = 6^d^ 


If  y^  is  chosen  so  that  y*^  =01'^+  E'^u^,  then  it  is  clear  that 
(y^,  u^;  TT,  6^)  is  a solution  for  (P),  (D).  Hence, 

Try*^  = S'^'d^  + TRn^' 

Repeating  this  argument  for  each  j £ we  get 
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(3.3) 


it(  E o^y^)  = 7r(  E o^d'^)  + Tr(  E aW)  ' 
j=l  ^ j=l  ^ j=l  ^ 

Again,  using  complementary  slackness  yields 

-_i-i  -i,i-i  -w  i-i 
-tB  z = V a z - z 

= v^a^  - A. (v.  + t. ) . 

1 1 1 

Summing  (3.3)  and  (3.^)  yields  the  result.  | 


O.h) 


The  proof  to  Lemma  3.3  shows  that  an  optimal  solution  to  (3.1) 
and  (3.2)  yields  profit  raaiximizing  production  vectors  by  letting 
y'^  = Thus  a proof  similar  to  that  for  Lemma  2.4  would 

demonstrate  the  following 


Lemma  3.4.  If  i £ m are  chosen  so  that  u^=0,  i£m  at  an 

optimal  solution  for  (3.1),  (3.2),  (tt,  z^  (i  G m),  y^  (d  € ^))  is  a 
competitive  equilibrium,  where  y^  = + Eu^. 

Also,  a proof  similar  to  Lemma  2.5  would  prove  the  key  property 


q . . . . q ... 

. = Tr(w^  + E c^cn^)  + v^a^  + E 6^(a^d^)  - A.  (v.  + t.)  > 0 
" d=l  ^ J=1  d ^ 1 X 


when  A, 't,  = 0 if  v.  < v*  where 
i i 11 


iiii^^  .i  i ^ i i. 

- <a,Bz  <w  + E(on)  + E‘^u‘^), 


V*  = max{y  z |z  > 0,  A z 


d=i 


D^u^  < o^d^,  for  all  j € £} 
J 


(3.5) 
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Next  we  show  how  the  freedom  to  add  constraints  to  the  households' 


L 


sets  allows  the  consideration  of  piecewise  linear  utility  functions. 


Suppose  u. (z^)  is  a concave,  piecewise  linear  utility  function  defined 
S.  ^ S. 


on  • Then  the  epograph  of  u^,  = { (z^,  z)  | z^  < Uj^(z),  z € 

is  a convex  set  which  has  a piecewise  linear  boundary.  Hence,  the 


epograph  of  u^^  can  be  written  as  a polyhedral  convex  set 


S. 


where  £ E ^ and  is  the  number  of  pieces  of  linearity  for 


u^.  Thus,  the  following  problems  are  equivalent 


maximum 


(5.6) 


subject  to  7r(B^z^)  < ttw^,  z^  > 0 


maximum 


subject  to  Zq  - g^jZ  < 


j 6 £, 


7r(B^z^)  < W-",  z-^  > 0 . 


i 1 


In  this  case  = (l>0, ...,0)  and  the  matrix  [A^|a^]  in  the  airxiliary 
linear  program  (3.1)  can  be  adjoined  to  the  matrix 

[e,  - g^^,  ...  , . 


In  the  next  two  chapters  we  describe  algorithms  which  utilize 
the  structure  of  the  models  which  were  described  here. 
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CHAPTER  III 


THE  BILINEAR  COMPLEMENTARITY  PROBLEM 
AND  AN  ALGORITHM  FOR  COMPUTING  ECONOMIC  EQUILIBRIA 

This  chapter  presents  a generalization  of  the  problem  described  in 
(II. 2. 3)  and  (II.P.^O  which  is  called  the  bilinear  complementarity 
problem.  This  approach  to  computing  equilibria  was  introduced  by 
Wilson  [1976].  In  the  first  section  we  essentially  reproduce  the  results 
of  that  paper  except  that  we  use  a new  proof  based  on  the  concept  of  a 
subdivided  complex  introduced  in  Chapter  II.  The  rest  of  the  chapter 
describes  a bilinear  complementarity  algorithm  in  enough  detail  to  be 
implemented  on  a computer.  The  results  of  computational  experiments 
using  this  first  implementation  of  a bilinear  complementarity  algorithm 
are  reported  in  Chapter  V. 


III.l.  The  Bilinear  Complementarity  Problem 

The  bilinear  complementarity  problem  (BCP)  is  to  find  x,  y > 0 
such  that 

Ax  + y = b,  A e 3R  b G ]R  ", 

(x,C^^)  _ x^y^  = 0 , i € m 

x^y^  = 0 , i e n \ m . 

Notice  that  if  m = 0,  the  problem  is  the  same  as  a linear 
complementarity  problem  (LCP)  (cf.  C.  Lemke  [I963]  or  B.  C.  Eaves  [1971b]). 
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It  is  easy  to  verify  that  the  equilihrixMi  problem,  as  we  have 
formulated  it  in  (2.5),  (2.4)  can  be  stated  in  this  fom.  It  is  an  area 
of  interest  to  determine  what  other  problems  can  be  formulated  as  a 
BCP  other  than  the  equilibrium  problem  or  those  which  can  be  formulated 
as  linear  complementary  problems.  We  present  the  general  problem  mainly 
to  ease  the  notational  burden. 

The  constraints  (4.1)  are  weakened  for  algorithmic  purposes  by 
adding  slack  variables  u^  > 0,  i S m to  the  bilinear  constraints; 

(x, C .)  - x.y.  = u. , i e m,  where  C £ 

Now  we  define  the  set 

W = {(x,y,u)  > 0|Ax  + y = b,(x,C  ^)+  x^y^  = u^,  i £ m,  x^,y^  =0,  i>m], 

(1.2) 

and  for  k = 0,1, ...,m  let  W,  be  the  subset  of  W for  which  u.  = 0 

for  i < k,  and  x^y^  = 0 for  i > k.  A solution  to  the  LCP  defined 

T 

by  (A,b)  yields  a point  in  if  x C > 0.  The  solution  to  the 

BCP  is  a point  in  W^.  The  bilinear  complementarity  algorithm  (BCA), 
to  be  defined  in  the  next  section  provides  a procedure  for  tracing  a 
path  from  Wq  to  provided  that  certain  assumptions  are  satisfied. 

We  shall  define  the  "algorithm"  by  defining  a mapping  F from 
an  n+m+l-dimensional  subdivided  complex  into  E and  proving  that 
F‘^(0)  is  a path  with  the  desired  properties  when  0 is  a good  value. 
First  we  define  a collection  of  sets  and  show  that  it  is  a subdivided 
complex. 
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A basis  P is  a subset  of  the  variables  i € ^ U^,  i £ m] 

(capital  can  be  considered  as  an  index  or  a name  for  x^,  etc.) 

consisting  of  (X^,Y^)  for  i < it  for  some  0 < k < m,  (Xj^,U^)  or 
(Y.,U.)  for  k < i < m,  and  exactly  one  member  of  (X^,Y^)  for  i > m. 

This  definition  is  slightly  more  restrictive  than  that  in  Wilson  [1976]. 

A path-basis  P(i)  for  i = k+1  is  obtained  by  adding  the  (k+1) 
element  to  a basis  3.  A sub-basis  3(i,j)  for  i ^ j is  obtained  by 
deleting  a j member  of  P(i). 

Now  if  3 is  a basis,  path  basis  or  subbasis,  then  ^(3)  is 
the  non-negative  orthant  of  indexed  by  the  set  of  variables  3. 

Note  th8,t  if  3 is  a basis  or  sub-basis,  then  (9'(3)  is  an  n+m-dimensional 
set,  and  if  3 is  a path  basis,  then  6?(3)  is  an  n+m+l-dimensional  set. 

/s 

Also,  a facet  of  3 is  an  orthant  corresponding  to  a basis  or  path 
basis.  Define  the  collection  of  cells 

s {<3'(3)|3  is  a path  basis]  . (l.3) 

The  followir.g  assumptions  are  of  vital  importance. 

Assumption  1.1:  In  W,  (x,C  ,)  > 0 for  each  i £ m. 

~ . , 1 

In  the  equilibrium  problem,  this  property  is  satisfied  by  Lemma 
2.3,  which  depended  upon  the  positivity  of  w^  and  the  choice  of  v^, 
i £ m. 
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Assumption  1.2;  The  subset  W*  = {(x,y,u)  G w|  Z u.  > 0)  is  bounded. 

i=l 


The  conditions  on  the  economic  problem  sufficient  for  assun^jtion 
1.2  to  hold  will  be  studied  in  the  next  chapter.  This  assiunption  implies 
that  we  can  confine  our  attention  to  the  compact  set 

K = {(x,y,u)  e < (x,y,u)  < ke) 

where  k>0  is  a sufficiently  large  constant. 

{^(P)  n k|^  is  a path  basis}  . (1.4) 

and  assume  that  there  are  no  redundancies  in  the  problem  so  that  each 
cell  C of  is  m+n+1  dimensional.  Define  M = Now  we  can 

prove  the 


Lemma  l.j.  is  a finite,  compact  subdivided  (nH-n+l) -complex. 

Proof.  The  number  of  path  bases  is  less  than  which  is  finite. 

is  closed  and  K is  compact  so,  by  (1.4)  each  cell  C £ 
is  a compact  (m+n+l)-cell.  All  that  remains  to  show  is  that  (M, ^) 
satisfies  the  definition  for  a subdivided  complex  (II.3.2,  Part  1). 

Since  there  are  only  a finite  n'omber  of  (m+n+l)-cells  in 
property  (c)  is  satisfied. 

By  the  definition  of  W,  the  facets  of  a cell  C of 
correspond  to  some  variable  x^  (or  y^  or  u^)  being  held  at  zero, 
where  £ 3 and  C = ^(3)  0 W*.  Thus,  since  we  can  associate 


2Q 


I 


each  cell 


Z with  a oasis  of  sub-basis,  we  can  associate 

each  cell  C c C c:  with  a set  of  variables  P not  con- 

strained to  be  zero.  Define  cp(C)  = P to  be  the  set  of  variables 
associated  with  C and  let  cp  ^(P)  = C be  the  cell  associated  with 

a)  (any  two  (m+n+1) -cells  of  are  disjoint  or  meet  in  a 

common  face).  Let  C^,  i = 1,2.  If 

3^0  3^=  , then  fl  = 9^,  a common  face.  Otherwise,  it  is  clear 

from  the  definition  of  face,  that  if  3^03^=  B,  then  C = cp  ^(3) 

is  a face  for  both  and  C^. 

b)  (Each  (m+n)-cell  of  ^ lies  in  at  most  two  (m+n+1) -cells. ) 
By  definition,  any  (m+n)-cell,  say  Cp  of  ^ is  the  face  of 

some  (m+n+1) -cell  of  , say  C^.  Suppose  3^  = 'P(C^)  and  3(k)  = 
cp(C^)  then  there  are  four  cases  to  consider  concerning  the  type  of 
variable  which  is  in  3(k)  but  not  32* 

1)  or  for  i < k.  In  this  case  there  is  no  path-basis 

which  can  be  formed  by  adding  a variable  to  3^.  Hence,  there  is  no 
(m+n+1) -cell,  other  tiiaji  C^,  which  contains  C^.  (Note  that  a SM 
in  this  case.) 

?)  Xj^,  Y^,  or 

2a)  If  X,  or  Y is  in  3(k)  but  not  3 , then  3„  is  a 

K K c 2 

basis.  By  definition  the  only  path  basis  which  can  be  formed  is  by 
adding  Uj^_^  to  3^5  let  3^  = 3^  U Then  = cp(3j)  is  the 

only  other  cell  containing  C^. 

2b)  If  (Uj^)  = 3(k)  \ 3j,  then  the  only  adjacent  cell  containing 

C^,  other  than  C^,  corresponds  to  the  path-basis  formed  by  adding  the 


50 


(k+1)®^  variable  ^k+1^  which  is  not  in 

3)  X^,  Y^,  or  for  k < i < m. 

5a)  If  or  is  in  6(k)  but  not  in  then  since 

G Pg,  the  adjacent  cell  is  identified  with  the  path-basis  formed 
tlr 

by  adding  the  i variable  which  was  not  in  P(k). 

3b)  If  is  in  p(k)  but  not  in  3^,  then  there  is  no 

adjacent  cell  to  which  contains  C^,  i.e.,  c:  dM. 

U)  X.  or  Y.  for  m < i < n. 

If  {X^}  = P(k)  \ then  where  ^ 

is  the  only  (m+n+l)-cell,  other  than  C^,  which  contains  C^.  The 

situation  is  analogous  when  &(k)  \ ^2  “ | 

Notice  that  the  adjacency  of  cells  and  indeed  the  definition  of 
(M,^)  is  dependent  upon  the  ordering  of  the  households  i = 1,2,...  ,m. 
Thus,  the  efficiency  of  the  algorithm  is  influenced  by  the  ordering  of 
the  households. 

Let  the  function  F = M K be  denoted  by 


F(x,y,u)  = 


Ax  + y - b 


(C  . ,x)  - x.y.  - u. , i £ m 

.,1  11  1 


F is  smooth  on  each  (m+n+l) -cell  C c 9/^7  . For  the  next  theorem,  we 
require  the 


Assumption  l.U.  The  cardinality  of  is  finite  and  odd. 
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Conditions  sufficient  to  ensure  that  Assumption  1.^4  hold  have 
been  supplied  by  various  authors  (e,g.,  B.  C.  Eaves  [1971b]).  Now  we 
prove  the  main  result. 

Theorem  1.9.  If  0 is  a good  value  for  F with  respect  to 

(M,  and  Assumptions  1.1,  1.2,  and  l.h  are  satisfied,  then 

i)  F ^(O)  is  a subdivided  1-complex  neat  in  (M, 

ii)  The  cardinality  of  the  set  of  solutions  to  the  BCP  is 
finite  and  odd. 

iii)  There  is  at  least  one  connected  component  of  F~^(0) 

with  one  boundary  noint  in  VL  and  the  other  in  W . 

U ra 

Proof,  i)  follows  immediately  from.  Theorem  II. 5. IT  (Part  1). 

Due  to  Proposition  11.5.18  (Part  1),  and  Assumption  1.1,  ii)  will  be 

demonstrated  if  it  can  be  shown  that  all  boundary  points  of  F ^(O) 

are  in  VL  and  W . 

0 m 

The  details  of  the  proof  of  Lemma  1.3  will  be  used  here. 

Remember  that  3(l)  refers  to  a path  basis  with  ^ 

t^’A^  = 0 for  i > 1.  If  the  face  of  = cp(p(0))  has 

^i'^i  “ then  by  case  2a),  is  the  only  cell  containing  C^. 

Hence,  c ciM.  The  union  of  all  such  facets  will  be  called  f*^. 

Similarly,  if  is  a face  of  = q)(B(ra))  such  that  u^  = 0,  i € m, 

then  is  the  only  cell  containing  (see  Case  2b).  These  facets 

in  the  boundary  of  M will  be  called  f®.  It  is  clear  from  the 
definitions  that  Wq  = F ^(o)  D f*^  and  ~ ^ ^(0)  f"”*  Since 


F ^(O)  is  neat  in  M,  if  we  can  show  that  F ^(O)  0 = VL  U W , 

0 m 

then  since  |Wq|  is  odd  and  |3f  ^(O) | is  even,  we  will  have  shown 

that  Iw  1 is  odd, 

' m 

From  the  proof  of  Lemma  l.i  it  is  evident  that  the  only  facets 
in  are  either  in  f^,  f"^,  or  of  the  type  described  in  1)  or  3b). 

Suppose  (x,y,u)  G F ^(O)  is  in  a cell  of  the  type  described  in  1). 

Then  x.  • y.  =0  and  0 = u.  = (x,C.  ) which  contradicts  Assumption 

1.1.  Suppose  (x,y,u)  € F ^(O)  is  in  a cell  of  the  type  described  in 
31)).  Again,  this  is  impossible  because  ® ' u^  = 0 is  a 

contradiction. 

Thus  F ^(O)  n Sm  = Wq  U and  ii)  is  proved. 

Since  each  connected  component  of  F ^(O)  with  a bc^j^idary  point 
in  f^  has  another  boundai-y  point  in  either  f^  or  f”',  a simple 
counting  argument  will  demonstrate  (iii).  I 


We  now  show  how  the  linear  pure  exchange  equilibrium  problem  can 
be  formulated  as  -a  BCP;  the  corresponding  formulation  for  the  model 
which  includes  production  is  analogous. 

The  linear  constraints  in  (2.3)  define  a system  of  n linear 
equations  in  2n  variables  Ax  + y = b where  n = m + i ^i  ^ ^ 

and  A is  a square  skew-symmetric  matrix  of  the  form 


i 

I 

I 


(1.5) 


and 
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The  variables  x,  y S have  the  form 


X = (i  e m),  tt;  (i  e m),  p] 

y = [t^  (i  e m),  Sj  (i  € m),  p] 

In  this  formulation  p is  constrained  to  be  non-negative,  so 

Tre  - p = 1 . 


in  the  dual  problem 


satisfied  until  a solution  to  the  BCP  has  been  attained. 

The  bilinear  constraints  (2.it)  are  expressed  by 

(C  .,x)  - x.y.  =0,  i € m 
.,1  11  — 

where 

C ^ = [0,...,  0,  0,  0 

th  ^ “ 

i element 

Corollary  1.6.  If  the  primal  and  dual  linear  programs  (V.2.1)  and 
(V,2.2)  have  unique  solutions,  then,  that  solution  (x^,y*^,  u*^)  £ Wq 
is  a boundary  point  of  some  path  y c F~^(0).  The  other  boundary 
point  (x*,  y*,  u*)  of  y is  a solution  to  the  equilibrium  problem 
(V.2.5,  2.4). 

Proof.  After  observing  that  any  point  in  Wq  is  a solution  to  the 
linear  programs  (2.1)  and  (2.2),  iii)  of  Theorem  1.5  shows  that 
(x*,y*,u*)  € W^.  Hence  ux-  = 0 and  the  constraints  (2.3)  and  (2.4) 
are  satisfied.  I 


III. 2.  The  Bilinear  Complementaritj 


jorithm 


In  the  remainder  of  this  chapter  we  explain  in  detail  an  algorithm 


for  following  the  path,  described  in  Section  1,  which  leads  to  an 
equilibrium. 
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The  algorithm  described  here  solves  an  equilibrium  problem  as 
formulated  in  Sections  2 and  3 of  Chapter  V,  not  the  bilinear  comple- 

1. 

i mentarity  problem  in  its  most  general  form.  The  reasons  for  this 

i. 

* restriction  are  twofold;  1)  we  have  not  found  any  other  classes  of 

I problems  which  can  be  formulated  as  BCP's,  2)  The  dimension  of  the 

linear  systems  of  equations  can  be  reduced  significantly  by  dealing 
with  D rather  that  A in  equation  (l.5).  We  still  refer  to  this 
algorithm  as  the  Bilinear  Complementarity  Algorithm  (BCA)  because  the 
complementary  property  ( II. 2. 7)  involving  and  u^  is  \rfiat 

motivates  the  definition  of  the  algorithm. 

In  the  remainder  of  this  section  we  give  a more  general  descrip- 
tion of  a bilinear  complementarity  algorithm.  In  the  following  sections 
we  describe  in  detail  the  linear  operations  we  use  to  reduce  the  dimen- 
sion of  the  path  following  subproblem  and  a specific  method  for  solving 
that  subproblem.  TVie  algorithm  we  describe  here  is  implemented  in  a 
computer  program  (BCA).  The  details  of  the  implementation  and  a report 
on  the  numerical  results  using  this  code  will  be  in  given  in  Chapter  V. 

Next  we  describe  in  more  detail  the  overall  structure  of  the 
algorithm.  This  structure  can  be  inferred  from  the  proof  that  (M,^) 
is  a subdivided  m-complex  and  the  proof  that  F '*'(0)  leads  through 
the  cells  of  ^ to  an  equilibrium  point  in  Section  1.  However,  to  make 
the  strategy  for  determining  which  cell  (or  path-basis)  is  adjacent  a 
little  clearer  we  will  now  give  a verbal  description  with  some  visual 
aids. 

A basic  point  in  W is  in  for  some  d = 0,1, ...,m.  It 

is  characterized  by  the  fact  that  u.  = 0 for  i < d and  u.  > 0 

1 — i 


•I 

I 


i = d+l,  . . . , n.  It  is  useful  to  depict  the  variables  which  are 
positive  and  zero  in  a chart: 


X 

y 

u 


123^^56789  10 
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+ 

+ 

0 

+ 

+ 

0 

0 

0 

+ 

+ 

a. 

0 

•r 

0 
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+ 

+ 

+ 

0 

0 

0 

FIGURE  2.1 


The  difficult  transitions  to  understand  are  those  for  which  the 
current  facet  which  f"^(0)  intersects  is  a basis.  Suppose  that  the 
last  cell  which  the  algorithin  passed  through  had  > 0 and  Figure  2,1 
represents  the  fact  that  u.^  .just  hit  zero.  The  next  step  is  to  intro- 
duce either  or  into  the  path-basis,  whichever  is  not  in  the 

current  basis.  In  this  case,  is  now  allowed  to  be  positive. 
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FIGURE  2.2 


Now  we  follow  F ^(O)  throiagh  the  cell  defined  by  the  current 
path  basis  until  some  other  variable  goes  to  zero.  We  hope  that  u^ 
hits  zero  so  that  we  move  towards  an  equilibrium  where  u^  = 0,  i £ m. 
Suppose  that  or  hits  zero  for  some  i > ^.  Then  the  adjacent 

path  basis  is  formed  by  adding  the  variable  complementary  to  the  one 
that  hit  zero,  i.e.,  if  yg  hits  zero,  introduce  Xg.  It  is  impossible 
for  x^  or  y^  to  hit  zero  for  i < 3 because  u^  = 0 implies 
x^*y^  > 0.  Similarly,  Uj^^  cannot  hit  zero  because  = 0 implies 

uj^  > 0 (recall  Il('’.7)).  There  is  one  more  possibility.  If  x^  or 
y^  hits  zero,  then  we  are  again  at  a basic  point,  and,  according  to 
the  definition  of  a path-basis,  the  only  adjacent  path-basis  is  formed 
by  introducing  Up.  Since  u^  is  no  longer  constrained  to  zero,  one 
could  say  that  we  have  "back-tracked"  in  our  goal  of  forcing  all  the 
budget  surpluses  to  zero.  However,  the  theory  of  Section  1 says  that 
this  goal  will  be  reached,  eventually.  In  fact,  by  the  proof  of 
Proposition  11.3.18  (Part  1),  only  a finite  number  of  cells  need  be 
traversed  before  an  equilibrium  point  is  reached. 

Below  is  a more  complete  specification  of  the  decisions  involved 
in  choosing  the  adjacent  cell. 

Bilinear  Complementarity  Algorithm 

0.  Initialize:  cUq  --  (x^,y^,u^)  is  determined  from  the  solution  of 

the  linear  program  ( II. 1) , with  the  associated  basis  m is 

the  number  of  consuners.  Let  j :=  0,  d :=  1. 

1.  If  d = m+1,  STOP  - coj  is  an  equilibrium  point.  Otherwise  let 

P.  :=  P.  U ((X  } or  (Y,)),  whichever  is  not  in  0..  Go  to  Step  L. 
J 0 u a j 
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2.  ii.  :=  P.  U U , go  to  4. 

J J 

3.  If  {Y  } = 3 \ 3 , then  3.  = 3 U {X  ] . 

i J J J ^ 

If  {X  } = 3.  A 3.,  then  3.  = 3.  U {Y  ) . 

J"-L  ' «J  ,J  J 1 

4.  Follow  the  path  (f|  )~^(0)  from  oo.  into  the  interior  of 

^,3.)  ^ 

J 

(^(3.)  until  the  cuiye  intersects  a facet,  call  it  <^(3..,). 

Let  CO.,,  he  the  point  in  that  intersection,  j :=  j+1. 

J+1 


5.  a)  If  CO.  e W,,  d ;=  d+1,  go  to  1;  else, 

J cl 

b)  If  co^  ^ 1^  :=d-l,  go  to  2;  else, 

c)  CO.  is  a sub-basic  point,  go  to  5. 

0 


III. 3.  Exploiting  the  LP  Structure 

We  shall  now  change  from  the  notation  of  the  bilinear  con^ile- 
mentarity  problem  to  a linear  programming  type  of  notation.  It  will 
be  important  to  distinguish  between  the  primal  and  dual  structural 
variables  contained  in  x and  the  primal  and  dual  slack  variables  in  y. 
The  piecewise  linear  equilibrium  problem  can  be  formulated  as 

follows: 

find  (x,t.  A,  O >0 


u,  = 


Dx  + It  = b 
AD  - /;i  = c 

(C  a)  - Ti^t^  =0,  i e m 

A^t^  = 0,  i € k\  m 

(x,0  = 0, 

5c> 


(4.1) 


such  that 


where 


X,  ^ G 


t,  A e ]R 


and  D G E 


kx  i 


It  is  clear  that  either  the  constraints  of  Il(r.ji)  and  (2.k) 
or  the  equilibrium  conditions  in  II. 5. may  be  expressed  as  above. 

The  size  of  k is  determined  by  the  number  of  consumers,  the  number  of 
goods  and  the  number  of  exogenous  constraints  on  the  consumers  and 
producers. 

The  fir.st  m components  of  the  primal  right-hand  side  are  the 
initial  utility  ‘'.evels  i G m.  To  determine  these  values,  one 

can  either  solve  m smaller  linear  programs  as  described  in  (2.6)  or 
(5.5)  of  Chapter  V or  let  v^  = 0,  i G m.  It  is  strongly  suggested 
that  the  former  course  be  taken  because  good  starting  values  will 
significantly  reduce  the  run  time  of  the  BCA,  and  the  choice  of  = 0 
will  usually  result  in  a highly  degenerate  solution  to  the  linear 
program,  which  may  cause  problems  with  the  BCA.  Once  these  parameters 
have  been  established,  the  initial  point  tiP  = is 

determined  by  solving  the  auxiliary  linear  program 


m&ximize  cx 

subject  to  Dx  + t = b (*<•?) 

X,  t > 0, 

to  find  (x*^,t^).  If  the  revised  simplex  method  (Dantzig  [I963])  is 
used,  the  optimal  shadow  prices  }P  and  relative  costs  can  be 

extracted.  Finally,  the  initial  budget  surpluses 


ho 


I 


i G m 


cem  be  determined.  By  the  complementary  olackness  condition  of  the 

optimal  primal  and  dual  variables,  A?-t^  - 0,  i C m,  and,  hence, 

0 0 

> 0 for  each  i G m.  Let  - be  the  index  set  of  the  k basic 


primal  variables  and  ^ ^ B \ ^ basic  dual  variables. 

Then  the  initial  (BCA)  basis  consists  of  (x,t)  the  primal  basic 

a 

variables,  {'h,0_Q>  dual  basic  variables,  and  u.,  i G m. 
c 

The  non-basic  variable  t^  will  be  added  to  iP  to  form  the 
-C) 

first  path  basis  p>  . The  path  which  the  BCA  follows  is  determined  by 


0 


(^l  ^ ) ^(O)  where  F(  is  defined  in  Section  1. 

(7(r) 

It  would  be  possible  to  follow  this  path  throxigh  the  (m+n+l)  dimensional 
cell  but  we  can  use  the  largely  linear  structure  of  (4.1)  and 

the  positivity  of  (c  .,k).  1 m,  to  reduce  the  dimension  of  the  path 
following  portion  of  the  algorithm.  In  a typical  path-basis  ? there 
are  n variables  which  are  basic  for  either  the  primal  or  dual  problem 
and  m+1  other  variables.  Tliere  are  d variables, where  d has  the  same 
value  as  in  the  description  of  the  BCA, among  (A^,  i G m),  (t^,  i G m) 
which  are  in  (P  but  not  in  the  linear  orogramming  basis,  call  them  A 

u 

Q 

and  t^.  Also  in  P are  the  m-d+1  bu'iget  surpluses  u^  which  are 

positive.  Using  the  primal  and  dual  basis  inverses,  we  can  represent 

the  n basic  variables  in  terms  of  the  d variables  A and  t , so 

U v’ 

we  have  reduced  the  dimension  by  n.  The  variables  u, u must 

d+i  m 

stay  positive  because  = 0,  i = d+1,  ..., 


Ul 


m,  so  we  can  ignore  them. 


The  variable  u,  and  the  constraint  = (C  ,,  a)  - A, ’t, 

d d . , d d d 

can  be  ignored  if  we  treat  (c  A)  - A^t^  >0  as 

a constraint  helping  to  define  the  cell  we  axe  in.  Thus,  we  have 
reduced  the  path  following  problem  to  one  in  the  d variables  A 
and  t and  the  d-1  equations,  (C  . , A)  - A.t.  = 0,  i € d-1, 

V .,111  

corresponding  to  u^  = 0,  i 6 d-1.  Next  we  describe  the  details  of  how 
this  reduction  of  dimension  is  carried  out  in  practice. 

In  linear  programming  codes  it  is  common  to  maintain  two  integer 
vectors  which  can  pick  out  the  basic  variables  and  tell  which  row  of 
the  matrix  they  pivot  on.  The  row  a basic  column  "pivots  on"  is  the 
index  corresponding  to  the  winner  of  the  min-ratio  test  performed  when 
that  col'imn  entered  the  basis.  We  shall  also  maintain  analogous  index 
vectors  for  the  dual  variables.  To  summarize,  the  following  information 
is  maintained  and  revised  during  the  course  of  the  BCA. 


<j(i;  = 
CT(i)  = 

r(J)  = 


-(a)  = 


the  primal  basic  variable  which  pivots  on  row  i,  i £ k 
the  dual  basic  variable  which  pivots  on  row  i,  i G i. 

( 0 if  a ^ CT 


a e n. 


I i if  a £ O'  pivots  on  row  i, 

I 0 if  j ^ 

I i if  a £ ^ column  a pivots  on  row  i, 


a G n. 


lj  = index  set  of  those  A-variables  not  in  u,  but  allowed  to  be  positive 
V = index  set  of  those  t-variables  not  in  cr,  but  allowed  to  be  positive 
d = the  number  of  pairs  (Aj^,t^)  which  are  both  allowed  to  be  positive. 


Remark,  a)  (j  U v = d 

b)  The  current  path  basis  3 can  be  written 

3 = {(x,t)  , (X.O  , A ,t  , u , u , ...  , u ). 

cr  - n V d d+1  m 

cr  ^ 

The  variables  X and  t will  be  referred  to  as  superbasic 

M V — * 

variables  (to  follow  the  terminology  of  Murtagh  and  Saunders  [1977] 
in  reference  to  the  nonbasic  but  positive  variables  in  their  GRG 
algorithm) . They  are  independent  variables  which  determine  the  values 
of  the  dependent  or  LP-basic  variables.  Next  we  show  how  the  dependence 
can  be  numerically  calculated. 

Let  the  primal  k x k basis  matrix  P be  partitioned  as 


X t 

CT  cr 


P = 


with  the  columns  of  P permuted  so  the  basic  slack  colcanns  are  on  the 
right  and  the  rows  they  pivot  on  are  at  the  bottom.  One  can  write 


5 = 


B 

E ■ 

P-"  = 

0 

A 

F 

and 

-ab"^ 

I 

(1^.2) 


Dividing  the  primal  system  into  basic  and  non-basic  variables,  we  have 


■ B 

0 

/ X \ 

■ E 

I 

1 -2.]  + 

A 

I _ 

F 

0 

If  we  want  to  see  the  effect  of  t^,  j G v on  the  basic  variables, 


^5 


ons  merely  updates  e^,  a column  of  I,  in  the  usual  manner.  That  is 
solve 

Ps  = e^  , Ph  = b , (4.3) 


so  that 


it)^ 


[S]  t + b 

vJ 


Clearly,  the  permutation  of  the  rows  and  columns  of  P was  not 
necessary  in  this  case,  because  we  are  assuming  that  the  primal  basis 
is  factored  so  that  it  is  easy  to  solve  systems  such  as  (4.5).  It  is 
useful,  however,  to  use  the  permuted  from  of  D to  describe  the  calcu- 
lations in  the  dual  system.  One  must  keep  in  mind  that,  in  practice, 

D is  not  actually  permuted;  the  index  sets  ct  and  v allow  one  to 
pick  out  the  elements  needed  to  perform  the  calculations  below.  We  write 
the  dual  system  in  terms  of  its  basic  variables  in  ct  and  its  non-basic 


variables  in 

cr. 

B 

E ■ 

A 

F 

(A_,^_) 

CT  cr 

_ 0 

-I  _ 

-I 

0 

= c = (c\c^) 


Multiplying  by  the  dual  basis  inverse. 


■ B 

E ^ 

-1 

0 

-I  J 

0 

-I 

yields  the  updated  tableau 
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+ (vy 

a CT 


ab‘ 

AB"'^-F 

-B-^ 

-b"^ 

= (c^"\  - c^)  . (U.U) 


Suppose  we  want  to  determine  the  effect  of  A.,  j £ upon  the 

0 

basic  variables;  i.e.,  we  want  to  calculate  the  appropriate  row  of 

(AB  AB  ^ - F),  Now,  -AB  ^ is  part  of  the  primal  basis.  Let  e. 

"th  T T 

be  the  j unit  vector  and  solve  s P - e..  Then  pick  out  the  components 
T 

of  s which  correspond  to  the  columns  of  P (a(i)  < £)  and 

call  the  result  s^.  To  get  the  appropriate  row  of  AB  ^-F,  just 
calculate 


-T  -T 

-S2  = - - F . . 


Then  we  have 


(A_,?_)  = + (c^B'\  c^B’^  - Cg)  , 


(4.5) 


If  we  want  to  update  the  row  corresponding  to  ^ . for  some 

j S a , the  procedure  is  very  similar.  We  want  to  calculate  the 

appropriate  row  of  (B  ^,B  ^).  The  desired  row  of  b"^  is  that  row 

corresponding  to  the  column  of  B associated  with  x.,  row  r(j)- 

J 

Suppose  this  is  the  q row  of  P.  Then  (B~  is  found  by  solving 

T„  T 
s P = e 

q 

T 

and  extracting  the  components  s^  such  that  a(i)  < i to  form 


-T  -T 

Let  s^  = s^E  and  we  have 


c^B’^  - Cg)  . 


(4.6) 


The  set  of  i^dated  columns  corresponding  to  t^  will  be 
k X I vl 

G1  £ ]R  ' ' and  the  transpose  of  the  updated  rows  corresponding 


to  A will  be  GC  £ 3R 
u 


The  following  equations  give  us  the 


basic  variables  in  terms  of  the  superbasics 


Git  + b > 0 

V — 


(4.7) 


= G?A  + c > 0 


We  have  described  how  G^  and  G^  can  be  computed  from 
scratchy  below  we  will  describe  how  G^  and  G^  are  updated  as  the 
basis  changes  and  superbasic  variables  are  added  and  dropped. 

First,  however,  we  describe  how  the  first  d-1  bilinear 
equations  (^.1)  are  expressed  as  f(A  ,t  ),  a function  of  the  super- 

U V 

basic  variables.  Only  the  bilinear  equations  which  are  binding 
(u^  = 0,  i £ d-1)  are  included  in  f.  Hence,  f maps  into 

The  inequality,  (C  ^,a)  - A^t^  > 0,  referred  to  on  page  42 
is  called  the  bilinear  inequality.  It  also  will  be  reduced  to  a 
function 


of  the  superbasics.  This  inequality,  along  with  the  inequalities  in 
(1.7)  determine  the  cell  which  the  algorithm  is  currently  concerned  with. 


h6 


By  substituting  the  correct  expression  in  (4.7)  for  the  basic 


primal  and  dual  variables  involved  in  the  first  d bilinear  functionals, 
vre  csui  write 

u = h^  + D A - diag(X  ) (F  t + e 1 
U■l-l^  plvl 

u = h-  + D-A  - diag(t  )(F  A + e^)  , ('+•9) 

y 2 2^  y2|i  2 


where  diag(y),  y € IR*^,  is  the  diagonal  nxn  matrix  with  y as  the  diagonal. 
Basically,  f(A  ,t  ) could  be  written  as  f(A  ,t  ) = (u  ,u  ),  but  we 

p V P V |J  V 

eliminate  the  functional  corresponding  to  u,  and  let  q(A  ,t  ) = u. 

CL  ^ V d 


as  it  is  expressed  in  (4.9). 

Now  we  are  in  a position  to  describe  a method  for  executing 
Step  4 of  the  BCA. 


III. 4.  The  Path-Following  Subroutine 

It  has  been  shown  that  the  path  defined  by 


f"^(0)  = {(A  ,t  )|f(A  ,t  ) = 0] 

p V p V 


in  the  cell  S 

J 


as  defined  in  (4.7)  and  (4.8)  is  identical  to  the  path 
described  in  Step  4 of  the  BCA.  Thus,  we  consider  an 


algorithm  for  following  f~^(0)  from  an  initial  point  (A^,t^)  £ 8S 

P V 

into  the  cell  S until  the  opposite  boundary  point  of  f ^(O)  D S 
is  reached. 

The  algorithm  to  be  described  here  is  an  adaptation  of  the  path 
following  algorithm  in  Section  III. 3.  The  adaptations  are  intended  to 


'*7 


speed  the  progress  across  a cell  on  the  assmption  that  a)  f~^(0)  is  nearly 

linear  and  b)  the  cells  S are  in  general  so  small  that  the  non- 

linearities  of  f ^(O)  are  not  significant.  The  linear  approximation 
“1  C 0 

to  f (O)  at  (A  ,t  ) is  used  to  make  a guess  at  which  facet  t of 

u V 

S f ^(O)  will  intersect.  Newton's  method  will  be  used  to  calculate 
the  intersection  of  x with  f ^(O). 

For  simplicity  let  the  variables  (X  ,t  ) be  represented  by 

u V 

z E E . Since  the  superbasic  variables  are  changing  with  every 
change  of  cell,  the  association  of  (A  ,t  ) with  z is  only  temporary.  In 
Chapter  III  (Part  1)  we  discuss  in  detail  how  the  tangent  to  f“^(0)  can  be 
calculated.  Suppose  that  Zq  £ SS  is  our  initial  point.  There 
calculate  which  satisfies 


f’(Zo)io  -0  , 


iUo''  ==  ^ ^ 


and  Zq  points  into  S. 

We  only  require  that  z^  have  norm  one  so  that  we  can  easily 
measure  distance  along 


T(a)  = Zq  + ^ a > 0 j 

the  linear  approximation  to  f ^(O)  at  Zq. 

The  next  step  is  to  determine  how  far  one  can  move  along 
T(a)  until  some  facet  t of  S is  hit.  For  the  linearly  defined 
facets  this  corresponds  to  a ratio  test  in  linear  programming.  We 
want  the  smallest  positive  a such  that 


t8 


i £ k , 


(Gl-T(a))^  + b.  --  0 , 

(GP-T(ar:).  + j i . 

Let  a*-  be  the  smallest  positive  root  of  the  n equations  above. 

There  is  one  facet  defined  by  the  bilinear  equation  q(z)  = 0. 

To  find  the  point  where  T(q)  first  intersects  this  facet  we  solve 
for  the  smallest  positive  root  of  the  quadratic  equation 

q(T(a))  = 0 . 

Let  a be  that  root,  and  replace  with  the  minimum  of  and  a. 

If  Q*  is  larger  than  some  maocimum  step  size  a , we  let  0^/2  be 
the  step  length  and  return  to  the  curve  f ^(O)  along  a normal  hyper- 
plane to  T(a)  as  in  Algorithm  III,3-1  (Part  1} . If  a*  is  less  than  a 

we  include  the  equation  defining  the  facet  containing  T(q*)  = z^ 

cl  d.  1 

with  the  d-1  functionals  in  f » ]R  “ -4  IR  and  use  Newton's  method 
to  solve  that  system  of  equations.  If  the  resulting  solution  is  in 
S (or  nearly  so),  the  desired  endpoint  of  f ^(O)  and  the  desired 
facet  of  S have  been  found.  If  z*  ^ S,  then  a point  z'  on 
the  segment  conv[ z^, z*]  is  found  which  is  on  one  of  the  violated 
constraints.  Newton's  method  is  a^ain  initiated  at  z^  to  find  the 
intersection  of  f ^(O)  with  this  facet  of  S.  In  practice,  it  is 
rare  for  T(a)  to  pick  out  the  wrong  facet,  and,  if  that  happens, 
the  procedure  described  above  usually  has  to  be  performed  only  once. 


U9 


steps ize  is  reduced 


1.  o*  > a , 

i max' 

2.  zi;  ^ S 

3.  z'  € convfz^jZ*]  is  determined. 

FIGURE  4.1 

Above  is  an  example  of  the  corrective  mechanisms  of  the 
algorithm  for  following  f”^(0).  To  find  the  correct  endpoint. 

It  required  2 tangent  calculations  and  8 Newton  iterations.  Our 
con^)utational  experience  suggests  that  for  most  applications  of  this 
algorithm  only  one  tangent  calculation  and  one  Newton  step  are  required. 


zorithm  4.2. 


0.  We  are  given  a,  4,  v,  d,  G^,  G^,  f and  q as  defined  in  Section  2. 

Let  z,  = (a  ,t  ) determine  the  initial  point  in  a facet  of  S. 

J-  1-1  V 

The  particular  facet  of  S is  determined  by  the  pair  (pd, s)  where 


pd 


If  f ^ -0  - 

if  q(Z]_)  = 0 , 


(4.1) 


Call  the  binding  constraint  b (z)  = 0.  , e^,  e,  > 0 are  fixed 

S i tf" 

parameters.  Let  i :=  1. 


1.  Calculate  f’(z^)  = [H|h]  £ . 

2.  Solve  Hy  = -h.  (Since  f'(z^  ) is  of  full  rank,  if  det  H = 0, 
choose  another  (d-1)  x (d-1)  submatrix  of  f'(z^),  H',  and  let 
h'  be  the  vector  left  over.  Set  H :=  H',  h :=  h'  and  repeat 
this  step.) 

3.  Let  CD  :=  sgn(det  H),  if  i >1  go  to  5. 

4.  6 :=  sgn  (Vb^(z^),  (y,l)), 

CO  :=  6'co,  go  to  6. 

5.  8 :=  CO-CO. 

6.  :=  5 , define  T^(a)  = z^  + z^. 
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w 

f 

i 

I 7.  Solve  the  n linear  equations  in  o.  defined  by 

<(Gi).^ T(a))  + = 0,  i e k 

<(G2)i^  , T(a))  + c.  = 0,  i € i , 

and  the  quadratic  defined  by 

q(Q(a))  = 0 . 

Let  be  the  minimum  positive  real  root  of  these  equations  and 

let  (pd,  r)  define  the  facet  which  Q(cd‘’)  is  contained  in  in  the 
same  manner  as  (4.1).  Let  b^(z)  =0  be  the  equation  defining 
' this  facet. 

Let  A = 1 and  ■P  = T(z*) 

8.  Find  the  minimum  positive  scalar  O.  such  that 

T^(a)  = 0 , i e d. 

If  a > q:^  go  to  9-  Otherwise,  let  r be  the  index  such  that 
T^(a)  =0,  and  let  pd  = 1 or  -1  depending  on  whether  z^ 

' corresponds  to  t on  A . Let  A = 2,  and  Cf**-  ;=  a,  z ~ T(o*), 

i ""  "" 

I and  b (z)  = z . 

• r r 

9.  If  < a go  to  Step  10. 

max 

9.5  Otherwise,  let  :=  or^/2, 

Z = T(q*), 

b^(z)  = (z^,  P-z)  , 

£ind  A = 5. 
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10,  Define 


until  II  < Gg.  If  the  termination  criterion  is  never 

realized  go  to  Step  9.5. 

^+1 

11.  If  A = 1 or  2 go  to  12,  else  A = 3.  Let  i :=  i+1,  ;=  z 

and  go  to  1. 

12.  Let  (z*)  = z"^^^  . Check  that 

((Gi)i^^,2*)  + b.  > , i C k, 

{ (G„) . ,z*>  + c > -e  , i e I , 

Cl  ±f  , 1 

and  q(z*)  > -c^.  If  z**-  satisfies  these  inequalities,  go  to  13. 
Otherwise,  define  pd  and  r to  correspond  to  the  most  infeasible 
constraint  and  represent  that  constraint  by  b^(z)  > 0.  Now  solve 

b {Oz*  + (l-G)z. ) = 0 

for  ft  £ [0.1].  Let  z^  = Qz*  + (l-9)z^,  A = 2,  and  go  to  Step  10. 

13.  If  pd  = 0 and  d = m go  to  Step  lU. 


‘=3 


Otherwise  STOP. 


lU.  To  get  a more  accurate  final  solution  calculate 


z*  ;=  z*  - g'(z*)"^  g(z*) 
aiid  repeat  until 

||g(z*)||  < 

and  then  STOP. 

We  shall  call  this  algorithm  the  endpoint  finding  subroutine. 

Several  remarks  can  be  made  about  this  subroutine. 

a)  The  theory  underlying  steps  1 through  6 is  contained  in  Section  III. 2 
on  the  orientation  of  paths. 

b)  The  principal  computational  effort  involved  with  this  subroutine 
can  be  divided  into  two  parts:  the  calculation  of  Jacobian  matrices 
and  the  solution  of  linear  systems,  and  those  operations  involving 
the  constraints  defining  S — the  calculation  of  or*-  and  checking 
that  z*  € S. 

c)  The  former  operations  are  of  order  (9'(d^)  while  the  latter  are 
(^(d.n).  If  n is  very  large  in  comparison  with  d (as  is 
usually  the  case),  then  the  latter  type  of  operation  taJtes  more 
time  than  the  former.  Since  d increases  as  the  BCA  runs  its 
course,  the  work  done  in  this  endpoint  subroutine  increases. 

d)  One  could  apply  Theorem  111.5.16  (Part  1)  to  state  that  if  is  chosen 

small  enough,  this  subroutine  will  follow  the  right  path  and 

converge  to  the  right  endpoint.  In  general  we  choose  a = 10  or  100 

nictx 


depending  on  the  problem.  The  Newton  iterations  cannot  be  guaranteed 
to  converge  with  such  large  step  lengths,  but  they  always  have. 


3)  Since  we  are  using  exact  partial  derivatives  in  Newton's  method  we 
can  state  that  these  subroutines  converge  quadratically  to  the 
endpoint  •j*  of  f ^(O)  (cf.  Ortega  and  Rheinboldt,  10.2.2, 
[1970]).  Tliis  means  that 


or  intuitively,  the  number  of  decimal  places  of  accuracy  eventually 
doubles  from  iteration  to  iteration. 

III. 5.  Moving  from  Cell  to  Cell 

An  important  factor  in  the  efficiency  of  the  BOA  is  the  matter 
of  how  quickly  all  of  the  revisions  can  be  made  to  the  basis,  updated 
columns,  and  function  definitions.  If  a good  linear  programming  routine 
is  used  to  solve  the  auxiliary  linear  program,  the  subroutines  from  that 
code  involved  in  updating  the  basis  and  solving  systems  will  perform 
the  cell  changing  operations  very  efficiently. 

A flowchart  of  the  decisions  and  operations  involved  in  deciding 
which  cell  is  adjacent  is  given  on  the  following  pages. 
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pd  :=  -1 

Choose  s,  the  incoming 
basic  variable;  Pivot  on 
and  G^;  Update  basis 

factorization  and  index  sets; 
(s  in,  r out); 

Superbasic  change: 

(\  in,  tg  out); 

Revise  definition  of  f; 

n • S T* 


from  super- 
basics, 
d ;=  d-1 
s : = d 
pd  : = 0 


go  to  1 


Remove  from  super- 
basics; Choose  s,  the 
incoming  basic  variable; 
Pivot  on  G^  and  G^. 
Update  current  basis‘s 
(s  in,  r out); 
Superbasic  change: 

(Ag  in,  tg  out) 

d :=  d-1 
s : = r 
pd  : = 0 


I 


yes 


\ 

[ 

I FIGUEE  5.1 


i 


Next  we  expeind  upon  some  of  the  operations  described  briefly 


in  the  large  boxes  of  the  flowchart. 

"Choose  s,  the  incoming  primal  (dual)  basic  variable." 

We  will  consider  the  case  when  pd  = -1.  What  transpired  during 

the  endpoint  subroutine  was  that  changes  in  A caused  the  constraint 

u 

G2  A + c > 0 
u * ~ 

to  become  binding,  where  r(3)  = ■^.  This  means  that  some  coefficient, 
G2..  =0,  i £ u.  Since  the  J variable  must  leave  the  dual  basis, 

AJ  X 

“til 

the  j ^ element  must  enter  the  primal  basis.  To  determine  s,  the 
variable  which  enters  the  dual  basis,  let 

s = arg  max | G2  | 

164 

We  choose  the  largest  element  in  absolute  value  to  aid  somewhat  in 
keeping  the  basis  matrix  well  conditioned,  because  G2.  will  be  the 

i/S 

pivot  element  in  a straight  forward  update  of  the  primal  basis  matrix. 

If  p = 1 and  £ = r(5)  we  similarly  find 

s = arg  max  1^1^^ | 
i £ V 

and  bring  the  column  corresponding  to  t^  into  the  basis  matrix  and 

til 

remove  the  column  corresponding  to  the  J ^ primal  variable. 

"Pivot  on  G1  and  G2" 


I 


Since  the  basis  is  to  change,  G1  and  G2  must  be  updated. 
Obviously,  one  could  update  the  basis  and  recalculate  G1  and  G2  as 
described  in  Section  1,  but  if  there  are  several  columns  in  G1  and  G2, 
this  would  be  rather  expensive.  Also  the  current  right-hand  sides 
b and  c must  be  updated.  We  will  describe  the  operations  performed 
if  p = 1 and  Z = r(r)  and  s = arg  max  |G1  . | has  been  chosen. 


First  we  pivot  to  update  Gl: 

Calculate  the  eta  vector  defined  by  Gl 


f.  :=  -Gl.  /G1„  , 

"1  is'  ’ 


i ^ i. 


Update  the  columns,  Gl  . and  S. 

• 0 


Gl^.  0 

Gl  . ;=  Gl  . + v-^  , 

•0  .J 

j 7^  s. 

(6.1) 

V :.S^ 

(6.2) 

b :=  b + V ^ . 

In  the  row  we  now  have 


E Gl  -t.  + 1-t  + b = 0 
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But  since  t > 0 
s - 


is  now  basic  we  want 


t 

s 


z 


Gl,.t. 
1-2  2 


+ b„  > 0 


so  we  must  change  the  signs  of  the  i 


th 


row's  coefficients: 


:=  ^ J s V,  j s (6.2) 

b • = -b  ^ . 

ijil 

Now  the  r dual  variable  is  entering  the  dual  basis,  but  since 

r > d,  this  is  not  a superbasic  variable.  So  we  must  calculate  the 

th 

updated  column  corresponding  to  the  r dual  variable  as  described  in 
(4.5)  or  (4.6)  depending  on  whether  the  variables  is  a A-  or  a 
^-variable.  Call  this  updated  col'jmn  Y. 

n = r(s)  is  the  row  we  are  pivoting  on  because  the  A^  is 
leaving  the  dual  basis.  Define  the  eta  vector 

’ 

‘ - 

as  usual  and  repeat  the  operations  in  (6.1),  (6.2)  and  (6.3)  with 
G2  replacing  Gl,  c replacing  b,  r replacing  s,  and  n replacing  £, 
"Update  current  basis  (r  in,  s out)" 


60 


I 
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Use  the  subroutines  of  the  linear  programming  code  to  update 
the  current  basis  factorization  when  column  j enters  the  basis  and 
column  s leaves.  Also  revise  the  index  sets  cr,  a,  y,  Y 'to  account 
for  the  change  as  follows: 

Let  k be  a dummy  variable  and 


k 

;=  r(s) 

k 

:=  f(r) 

c(k) 

;=  r 

a(k) 

: = s 

r(s) 

:=  0 

?(!•) 

:=  0 

r(r) 

:=  k 

r(s) 

:=  k . 

''Superbasis  change  (a  in,  t out)" 

s s 

Update  the  index  sets; 

p : = U U s , 
V :=  v\  s. 


Add  the  column 


corresponding  to  as  calculated  in  (4.5)  to  the  matrix  G2,  using 

the  new  basis  factorization  to  calculate  s^.  Remove  the  column  of  G1 
corresponding  to  t^ . 

"Revise  the  definition  of  f(A  ,t  )" 

I , V 


Actually  we  also  revise  the  definition  of  q(A  ,t  ),  the 

U V 

bilinear  inequality,  here  too.  This  is  just  a recalculation  of  the 
bilinear  equations  of  (1.9)  using  the  updated  matrices  G1  and  G2, 


6l 


( 

f 

and  the  revised  ri^t-hand  sides  b and  c to  again  express  the  basic 
variables  in  terms  of  the  new  superbasics. 

The  major  work  done  in  the  change  from  cell  to  cell  is  the  pivot 
on  G1  and  G2,  the  calculation  of  two  new  superbasic  columns,  and  the 
update  of  the  basis  matrix.  Clearly,  the  efficiency  of  this  portion 
of  the  algorithm  depends  upon  the  efficiency  of  the  particular  sub- 
routines of  the  LP  code  which  is  used. 

It  appears  to  be  that  an  interesting  area  for  algorithmic 
research  and  experimentation  lies  in  studying  the  endpoint  algorithm 
itself.  Many  variants  of  Newton' s method  could  be  used  to  solve  the 
nonlinear  equations  involved.  Perhaps  sophisticated  differential 
equation  methods  could  be  used  to  follow  the  path. 

In  Chapter  V we  shall  report  on  the  numerical  results  from 
the  testing  of  a code  which  implements  the  algorithm  described  in  this 
chapter. 
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CHAPTER  rv 


A HOMOTOFY  METHOD  FOR  COMiPllTIHG  EQUILIBRIA 
IV. 1.  A Convergence  Theorem 

The  bilinear  complementarity  algorithm  described  in  the  last 
chapter  is  essentially  a systematic  procedure  for  balancing  the  budgets 
of  the  m consumers  sequentially.  One  might  guess  that  a rather  large 
number  of  cells  would  be  traversed  before  equilibrium  is  reached.  For 
example,  at  least  m cells  must  be  traversed  in  which  the  curve 
F ^(O)  hits  a budget  constraint.  In  this  chapter,  we  investigate  a 
path  method  which  immediately  relaxes  "t^  =0,  i £ m,  and  adjusts 
the  variables  so  as  to  solve  for  all  of  the  budget  surpluses  at  once. 

We  use  the  homotopy  retraction  method  from  Section  IV. 1 (Part  1)  to  motivate  the 
construction  of  the  deformation  which  defines  the  path  of  interest. 

We  will  again  utilize  the  pure  exchange  model  of  Section  II. 2 as 
the  generic  example  to  ease  an  already  ciunbersome  notational  load. 

The  more  general  economy  of  Section  II. 3 can  be  dealt  with  in  an  analogous 
manner.  The  assumptions  which  guarantee  that  equilibria  exist  and 
coincide  with  quasi-equilibi ia,  II. 2. 3,  will  be  in  effect  here.  The 
utility  values  v^  vnll  be  chosen  in  accordance  wdth  (2,7).  The  same 
primal  and  dual  auxiliary  linear  programs  (II.2.1)  and  (2.2)  will  again  serve 
as  a vehicle  for  treating  the  equilibrium  problem  as  an  equation  solving 
problem. 


= (tt,  A,  ^^(ie  m),  p;s,t,z^(i  e m )>p) 


vAiich  satisfies 


where 


f (cp)  = 0 , 


CD  € D , 


(1.1) 


f(tD)  = 


tt(co)w1  - A^(cD)(t^(a))  + v^) 


'7r(a))w  - A (a))(t  (cd)  + v„) 

m tn  m 


(1.2) 


Tr(a>)  refers  to  the  first  m components  of  co,  etc.,  and 


D = {cd| 


r z - t.  = V.  , 
11 


Z + pe  + s = Z 

i=l  iem 

ttB^  - A^r^  - = 0 , 

Tre  - p = 1 
ITS  = 0 
= 0 , 
CD  > 0 , 


i € m 


i € m 


i 6 m 


(1.3) 


Since  f maps  (Tr,A,t)  € into  ]R  ”,  it  does  not  appear 

that  the  path  methods  such  as  the  homotopy  retraction  or  strong  path 
methods  can  be  used  to  solve  this  problem.  But,  under  some  quite  reason- 
able assumptions,  it  can  be  shown  that  D is  a collection  of  m-dimensional 
polyhedral  sets  which  form  a subdivided  m-complex.  By  a careful  choice 
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of  an  initial  point,  an  analogue  of  the  homotopy  retraction  method  will 
be  used  to  define  a path  leading  to  an  equilibrium  point. 


In  the  definition  of  D (1.3)  the  complementary  relation 
A^‘t^  = 0,  i £ m,  is  not  enforced.  For  this  reason,  none  of  these 
variables  will  be  constrained  to  zero  during  the  course  of  the  algorithm. 
By  V.2.7  we  have  that  if  or  t^  = 0 then  f^(cD)  >0,  i £ m.  Thus, 

the  retraction  function  which  maps  points  in  ]R^  into  can  be 

defined  in  terms  of  aCoj)  £ or  t(a))  £ In  either  case,  the 

boundary  condition  that  f(cQ)  points  into  E™  for  A(u))  £ Se'^  (or 
t(a))  £ SE^)  will  be  satisfied.  Since  D is  not  a bounded  convex  set, 
we  must  make  some  assumptions  to  allow  us  to  conclude  that  the  path 
defined  by  the  deformation  is  bounded.  The  assumptions  required  are 
more  natural  when  the  A-variables  are  kept  nonnegative  by  the  intrinsic 
properties  of  the  path,  rather  than  the  t-variables.  The  sign  of  the 
t-variables  will  be  unrestricted,  but  if  a solution  of  (l.l)  can  be 
found,  then  f^(a')  = 0 implies  that  A^(ai)  and  t^(fo)  are  positive 
by  2.7  of  Chapter  II. 

The  preceding  discussion  motivates  consideration  of  the  set 
K 3 D,  defined  below.  If  we  define 

0)^=5  (z^(a)),  i £ m,  p(co),  7r(cu)), 

00“"=  (^^(cu),  i £ m.  p(co),  5(03))  , 

and  let 
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0 


^IT  1 
-B  r 


1 2 

A set  L = (o)^,  oi^,  A^,  t^)  of  n variables  (n  = Es^  + & +m+l) 

1 2 

consisting  of  exactly  one  of  each  pair  (05^,05^),  i € n-m  or  (A^,t^), 
i € ^ is  a feasible  basis  if  the  col’omns  of  [Ajl]  corresponding  to 
E ([A|I].j,)  is  a nonsingular  matrix  and  the  system 


[A|I], 


+ [A  ] = b , 

.4  .V  t^ 


1 2 . - 

0)  I CO  >0 

a a - 


(1.5) 


has  a solution,  where  A = {A. |A.  ^ Z)  and  t s (t  |t  ^ E] . 


A characteristic  set  of  variables  is  a set  T = Z U ^ U v for  some  feasible 
basis  Z.  Note  that  a basis  determines  a characteristic  set  (c-set), 
but  not  vice  versa.  Since  a c-set  r always  contains  all  of  the  A- 


and  t-variables,  it  can  always  be  determined  by  r = r\  {A^,ti,  i e m) . 
Let  f be  those  variables  complementary  to  y. 

Corresponding  to  each  c-set  of  variables  T,  we  define  the  set 


C.,  = <(d|[A  I] 


cu 


05 


T J 


05 


05^  > 0, 


(1.6) 


where 


0 1 


o 


and 


1. 


are  the  coliomns  corresponding  to  the  A-  and  t-variables.  These 

sets  will  satisfy  the  definition  of  a cell  {11.3,  Part  1)  if  we  assume  that 

each  cell  corresponding  to  a characteristic  set  of  variables  satisfies 

the  constraint  qualification.  This  will  imply  that  there  is  a point 

05  C Cp  such  that  » 0 and  » 0.  Let  ^ = {C^|r  is  a c-set]. 

(We  will  occasionally  refer  to  elements  of  ^ as  characteristic  cells 

or  c-cells  to  emphasize  their  correspondence  with  some  c-set  of  variables.) 
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Define  K = U C 

Cpe^t  ^ 

An  important  step  in  applying  the  theory  of  Chapter  II  (Part  1) 
is  to  show  that  (K, ^ ) is  a subdivided  m-complex. 

Lemma  1.1.  (K, is  a subdivided  m-complex. 

Proof.  We  make  use  of  the  one-to-one  correspondence  between  c-cells 
and  c-sets.  Property  (c)  if  1.5.2  follows  because  there  are  only  a 
finite  number  of  characteristic  sets. 

The  proof  follows  the  lines  of  Lemma  111.4.5  (Part  1)  because,  again, 
any  cell  C € ^ can  be  associated  with  a set  of  variables  (t^p,cOp) 
not  constrained  to  be  zero.  Define  cp(c)  =3  to  be  the  set  of  vari- 
ables associated  with  C and  let  cp  ^(3)  = C be  the  cell  associated 
with  3.  Property  (a)  now  follows  from  Part  a)  of  the  proof  of 
Lemma  111.4.5.  Property  (b)  follows  from  Case  4)  of  part  b)  of  the 
proof  of  the  same  lemma.  | 

Below  we  will  expand  the  subdivided  m-complex  (K,  ^1^)  by 
adding  the  parameter  G needed  to  define  the  deformation.  Clearly 
if  y = {cr|a  = c^  X [0,1],  '){)  and  L = a,  then  (L,^) 

is  a subdivided  (mtl) -complex.  Define  F:L  as 

F(cd,0)  = 0f(a))  - (l-0)(A(cn)  - A°)  (1.8) 
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It  will  be  assumed  that  - 'h(uP)  for  some  aP  d K 
which  is  easy  to  compute.  For  theoretical  purposes  it  will  be  useful 
to  choose  fo*^  such  that  hioP)  = 0 d This  appears  to  be  necessary 

in  order  to  show  that  F ^(O)  is  boimded.  The  process  of  following 
F~^(0)  to  an  equilibrium  point  will  be  called  the  hornotopy  retraction 
method. 

Next  we  describe  the  assuniptions  necessary  to  prove  that  we  can 
consider  a compact  subset  of  L. 

Assumption  1.2.  (Bounded  utility).  Even  household's  utility  is  bounded 
for  a fixed  level  of  exports.  That  is,  for  fixed  p and  any  i 6:  ra, 
the  linear  program 

. . i i 

maximize  y z 

subject  to  E B^z^  + pe  < Z wl  z^  > 0, 

ifm  i£m  ^ 


has  a finite  objective  value.  V/e  also  make  the  assumption  that  the 

utility  is  bounded  below  for  fixed  p. 

The  first  part  of  the  assumption  is  fairly  standard  in  the 

literature  (cf.  Dantzig,  Saves,  Gale  [1976]).  The  second  part  seems 

as  innocuous  as  the  first.  Even  though  there  might  be  an  unattractive 

activity  j C S.  for  some  i 0)  it  seems  reasonable  to  assume 

^ J 

that  household  i can  perform,  only  a finite  amount  of  that  activity, 
even  if  the  resources  of  the  other  households  are  made  available  to  him. 


j 


I; 

f 

t 


i 
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Asstunption  1.:^.  (Bounded  production) 

a)  The  amount  of  exports  the  system  is  capable  of  producing  is  boundedi 
the  optimal  value  of 

maximize  p 

subject  to  Tj  B^z^  + pe  < E z^  > 0 . 

i£m  iEm 

is  finite  (call  this  optimal  value  p^) . 

b)  The  solution  to  this  linear  program  and  its  dual  is  unique. 

Part  b)  will  not  necessarily  be  satisfied  for  an  arbitrary 
problem,  but  by  perturbing  the  objective  vector  and  the  sum  of  endowments 
this  property  can  be  satisfied. 

Assumption  l.U.  (7r(co),  w^)  >0,  Vi  £ m for  o)  G K. 

This  assumption  is  a weaJcening  of  Assumption  IV. 2. 5 (Part  1),  but  it 
also  implies  the  existence  of  equilibria  and  that  equilibria  and  quasi- 
equilibria coincide  (cf.  Arrow  and  Hahn  [1970]).  In  models  with 
production,  this  assumption  can  be  weakened  to 

, w^  + Z cr^cu^)  >0  V (.0  £ K. 

d=l  ^ 

This  assumption  implies  that  the  value  of  each  consuimer's  assets  is 
positive  for  any  realization  of  to  £ K.  This  condition  is  very  difficult 
to  verify  for  any  particular  problem,  in  contrast  to  the  assumption 
that  w^  » 0,  for  any  i £ m.  However,  computational  experience 


70 


t 

L 

f 

i 

I 


'I 


I 

E 


L 


indicates  that  Assamption  l.ii  is  often  satisfied  in  practice. 

Now  we  are  in  a position  to  prove  that  major  result  of  this 
section,  that  the  component  of  F ^(O)  leading  from  {aP ,0)  terminates 
at  a point  (m*, 1)  which  is  an  equilibrium  point.  First,  we  must  show 
how  to  choose  an  initial  point  lP  such  that  A(a)*^)  0 and 

{uP,0)  e 1. 

Let  (zq,  i e m,  p^,  s^)  solve  the  problem 


maximize 


subject  to  E + pe  + s = E w^ 

iGm  i€m 


(1.9) 


p,  s > 0,  z > 0,  iGm 


and  let  iGm,  p^,  tt^)  solve  the  dual  problem 


minimize 


subject  to 


7r(  E w^) 
iGm 

Tre  - p = 1 
= 0 , 

O,  TT  > 0,  > 0 


iGm 

iGm 


1 1 


m 


Then  let  t^^^  = y z^  - v^,  iGm  and  Aq  = 0 G IR  . It  is  easy  to 
see  that  (co°,0)  = (tt^,  A^,  (i  G m),  0^;  s^,  t^,  z^  (i  G m),  p^jO) 

is  the  only  point  in  L with  A(oo)  = 0 due  to  assumption  1.3(b). 

Our  deformation  F:L  -»1r”  has  the  simple  form 


i 
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F(a),e)  = efCoj)  - (i-e)  Mo^)  , 


(cu,  0)  e L . 


We  will  be  assuming  that  0 is  a good  value  for  F w, r,t.  (L,  ^ ). 


Lemma  1.6.  If  P is  the  component  of  F"^(0)  containing  aP , and 
Assvunptions  1.2  and  1.3  are  satisfied,  then  P is  a bounded  set. 


Proof.  Suppose  that  P(a);[0,T]  -» L is  a continuous  parametrization 
of  P where  T < + » . First  we  show  that  >v(P(a) ) » 0 for  any 
a > 0.  Since  we  begin  at  the  facet  of  L for  which  0=0,  it  is 
clear  that 


a=0 


> 0 . 


Since  A(0)  = 0,  Assumption  l.U  gives  us  that  f(co*^)  » 0. 
Also,  we  can  express  f solely  as  a function  of  A(u3),  initially, 
because  the  t variables  are  basic  at  the  start.  Thus, 


F*(a)°)  = [df'(A)  - (l-0)l|f(A)  + A]  . 


At  a = 0,  given  that  0(a)  > 0,  it  must  be  true  that  A(0)  satisfies 


-IA(0)  = -f(A(0)) 


or 


dA(a) 

da 


a=0 


= f(A(0))  » 0 . 


Hence,  there  is  some  e > 0 such  that  for  any  a £ (0,e],  A(a)  » 0. 
Suppose  there  is  an  a £ (€,T]  such  that  Aj^(a)  < 0 for  some  i £ m. 
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( 

i 

Then  by  the  continuity  of  A^(a),  there  must  be  some  a e (g , a]  such 
that  A^(q:)  = 0.  Then,  since  P(a)  € F~^(0), 

e(a)  •f^(a)(a) ) = o . 

This  is  impossible  because  (a)  9(a)  = 0 implies  that  A(a)  0 and, 
hence  P(0)  = P(a)  which  contradicts  the  fact  that  0 is  a good  value 
for  F,  and  (b)  A^(5)  =0  implies  f^(a)(S))  >0.  Thus  A(a)  »0 
for  any  a > 0. 

The  fact  that  A(qi)  » 0 for  any  a £ (0,T]  immediately  yields 
that  f(co)  > 0 for  any  oo  £ P because 

^i^“^  " ^ ’ i £ m . 

From  equation  (IV.2.5,  Part  1),  we  have 

m 

p(cd)  = Z f . (co)  > 0 , 00  £ P . 

i=l 

Thus,  by  Assumption  1.3(a),  0 < p(oo)  < p’^  for  any  go  £ P. 

It  is  well  known  that  the  optimal  value  of  a linear  program  is 
a continuous  function  of  the  right-hand  sides  for  which  the  linear 
program  has  a feasible  solution  and  finite  optimal  value.  Consider 
the  following  continuous  functions  of  p. 
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i £ m 


(J  (p)  = inax{r  z|ZBz<Ew-  pe,  z^  > 0,  i € m), 
i=l  i=l 


. . m . ® i 

cp  (p)  = min{r^z^|  Z B^z  < Z w - pe,  z^  > 0,  i £ m],  i £ m. 
i=l  i=l 


By  Assumption  1.2,  $j^(p)  and  cp^(p)  are  finite  for  any  p £ [0,p^],  and, 
by  continuity,  $^(p)  attains  its  maximum,  t^,  and  cp^(p)  attains  its 
minimum,  t^,  on  [0,p^]  for  any  i £ m.  From  the  definition  of  D ri  P 
(1.3)  we  have  for  oj  £ P. 

t7  - V.  < t.  (co)  <t'!'-v.,  i£m. 

1 1—1—1  i'  — 


Next  we  show  that  A(a))  is  bounded  for  cu  £ P. 

7r(cja)  £ {ttItt  >0,  -rre  = 1),  a compact  set.  Clearly  bj^(7r)  = ttB^ 

is  a continuous  function  of  v.  Thus  b^  is  bounded  above  by  k'e^ 

for  some  k > 0 and  e £ ]R  ^ is  a vector  of  ones. 

Since  is  an  insatiable  utility  function,  > 0 for  some 

0 

j £ s,.  Pick  J = arg  min('irf|rv  > 0)*  Then  the  constraint  A.y^  < ttB^ 

X 1 ^ & X i “ 

— k £ s. 

1 

< ke^  implies  that  < k/y^  = A^  < » for  any  i £ m.  Since 
A^  > 0,  {A(cd)|iio  € P)  is  a boiuided  set.  Since  t(cD)  and  ACcu)  are 
bounded  in  their  respective  subspaces,  we  can  use  (1.5)  to  claim  that 
all  other  variables  are  bounded.  Given  a feasible  basis  E,  we  have 


1 

00 

o 


2 

00 

0 


t 

a 


b - [A  |I  ] 

• U -v 


7^+ 


1 2 

Since  m and  oj  are  continuous  functions  of  A(m)  and  t(cD)  for 
a a 

each  of  the  finite  number  of  bases  z,  we  conclude  that  is  bounded 

for  any  i G 2n  when  oi  G P.  Therefore  P is  bounded,  | 

Once  again  we  redefine  our  subdivided  (irH-L) -complex.  Choose  a 
constant  Q > 0 large  enough  that  for  any  od  G P,  G (-Q,Q)  for  any 
i G Define  C s C H { (cn,  d)  |cd^  G [-Q,Q],  0 G [0,1]]  for  any  cell 

C E aC . Let  ^ be  the  collection  of  all  such  cells,  and  let 

M = U C.  Then  (M,  is  clearly  a subdivided  (m+l) -complex  which 

CG7)7  ' 

is  compact.  We  can  now  state  and  prove  the  main  result. 

Theorem  1.7.  If  0 is  a good  value  of  P with  respect  to  (M, 
and  Assumptions  1.2,  1.5,  and  1.4  are  satisfied,  then  the  component  of 
f”^(0)  containing  (a)*^,0)  leads  to  the  boundary  point  (a>*,  1)  and 
CO*  is  an  equilibrium  point  for  the  piecewise  linear  economy. 

Proof.  The  facets  of  the  cells  of  which  are  contained  in  the  boundary 

of  M are  those  corresponding  to  co^  = + Q,  i G 2n  or  the  0=0  and 

0-1  boundaries.  By  Corolla,ry  II. 5. 19  (Part  1),  the  opposite  boundary  point  to 

(to^,0)  in  P must  be  in  the  boundary  of  M.  Here  we  are  using  the 

fact  that  0 is  a good  value  for  the  compact  complex  M.  We  have 

already  shown  that  0(a)  > 0 for  a > 0 and  the  path  P never  hits 

a boundary  corresponding  to  (cu^j  = Q for  any  i G Thus,  P must 

lead  to  a boxindary  point  (cu*,l).  By  the  definition  of  F,  f(co*)  = 0. 

Also  the  proof  of  Lemma  1.6  shows  that  A(a)*)  » 0,  so  by  Assumption  1.4, 
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f 


1 


I 

i 


t((J0*)  » 0.  Finally,  since  cd*  G K,  we  appeal  to  Lemma  IV. 2. it  to 
conclude  that  cd*  is  an  equilibritjm  point.  I 


IV. 2.  The  Homotopy  Retraction  Algorithm  For  Solving  the  Piecewise 
Linear  Equilibrium  Problem. 

In  these  sections  we  describe  how  the  theory  of  Section  1 can 
be  implemented  as  an  algorithm  with  some  nice  convergence  properties. 
Much  of  the  description  here  depends  upon  methods  and  notation  that  was 
introduced  in  the  previous  section  of  this  chapter. 

We  are  putting  the  problem  into  the  form 
find  (x,t;A,  0>0 

such  that 


Dx  + It  = b 
AD  - U = c 


A.  -t.  =0  , 
1 1 


i € k\  m 

U,0  = 0 . 


k X k X m 

Where  again  D € ]R  , and  C G E contains  the  appro- 

priate constants  determining  the  income  of  the  consumers  as  discussed 
in  Lemma  III. 5. 5 (Part  1). 

Again,  we  reduce  the  problem  to  deal  only  with  the  superbasic 
variables  A and  t . We  will  write  the  defining  inequalities  for 

b V 

a typical  cell  in  terms  of  the  matrices  G1  and  G2  and  the  updated  right- 
hand  sides  calculated  as  in  Section  Thus  ^+>7  ) of  Chapter  III  can  be 


j 


written 
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The  initial  point  03°  = is  determined  by  solving 

the  auxiliary  linear  program  II. 3.1  rather  than  the  procedure  suggested 
in  Chapter  III  which  results  in  an  initial  point  cn'  such  that  ) =0, 

i £ m.  Solving  the  auxiliary  linear  program  will  result  in  a point 
which  is  much  closer  to  the  final  solution  than  a solution  for  which 
Xj^(u3)  =0,  i £ m.  Thus,  for  reasons  of  algorithmic  efficiency  we  choose 
an  initial  point  which  is  not  theoretically  guaranteed  to  converge. 
However,  in  nearly  all  examples  which  have  been  run  on  the  computer 
convergence  was  achieved. 

The  method  still  makes  use  of  some  of  the  desirable  feasibility 
properties  of  the  homotopy  retraction  method  by  calculating  the  boundary 
point  which  f(a5*^)  points  at  from  i £ m.  To  do  this  we  let 

a = min  {'K AuP)  /t.{oP))  , 
i£ra 

and  - ctf^(a))),  i £ m.  Then  £ 3Ir”  is  the  "initial 

point"  of  the  algorithm  even  though  the  algorithm  never  was  at  a point  oi 
such  that  ^(03)  = 

The  deformation  which,  along  with  the  subdivided  m-complex 
K.^),  defines  the  path  the  algorithm  follows  is  given  below. 
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F(a),0)  = 0f((D)  - (l-0)(A(co)  - A°) 


0°  = a/(l  + a) 


03  e K,  0 e [0,1]. 


Again,  the  variables  in  03  are  expressed  in  terms  of  A^,  and 
as  in  (1.9). 


F(o3,e)  = 


e(h, 


+ D,  A 


diag(A  )(F  -t 

(0  i V 


+ e ) - (1-0)  (A  - a"") 
j-  HP 


e(h  + A - diag(t  )(F  A +e) 

d d ^ A d ^ d 

- (1-0) (G2  A + c - A°) 

p.  U M p 


(5.1) 


It  is  important  to  choose  the  parameters  v*,  i £ m be  chosen 
large  enough  that  A^(to^)  > 0.  Often  the  choice  of  suggested,  in 

Section  III. 3 will  result  in  positive  dual  multipliers,  but  if  not, 
v^  must  be  increased  (parametrically)  until  all  A^  are  positive, 
i = 1,2,...,  m.  If  some  A^(o3®)  = 0,  then  a = 0,  6*^  = 0,  and 
A^  = A(o3^).  In  the  optimal  solution  to  the  auxiliary  linear  program, 
the  slack  variables  on  the  first  m rows,  t^,  i 6 m are  never  basic, 
so  t^  = (t^,t2,  ...  , t^),  and  p = (Z.  Thus,  for  the  first  cell,  the 
function  is 


and 


F(a',0)  = 0(h  - diag(t  )(e^))  - (1-0)  (5  - A°) 

^ V 2 pi 


F’  ((D,  0) 


0=0 


[-0  diag(e2)|h2  - diag(t^)e2]  . 
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1 


But  if  9^  = 0,  the  matrix  F'  (cd,  6)  „ is  clearly  not  of  full  rank, 

and  hence  the  first  step  of  the  path  following  algorithm  cannot  be 
executed.  Thus,  the  importance  of  choosing  so  that  A^(fJD^)  > 0, 

i e m is  clear. 

In  Section  II. 1 (Part  1)  on  the  discussion  of  the  homotopy  retraction 
method  the  dimension  of  the  problem  was  increased  by  one  by  the  homotopy 
parameter  0.  This  was  done  primarily  for  theoretical  reasons — it  is 
easier  to  analyze  the  behavior  of  a path  \diich  ends  at  a boundary  point 
of  a complex,  rather  than  a point  at  which  the  deformation  function  is 
undefined.  However,  one  might  suspect  that  it  would  be  worthwhile  to 
formulate  the  deformation  in  the  lower  dimension  for  computational 
purposes.  For  example,  if  = 0,  the  deformation  could  be  defined 

A (m) 

h(a>)  = f(cD) 

and  the  path  would  be  h However,  f^(a))  ->  0 and  X^(cn)/f^(cD)  ->  + oo 

as  the  algorithm  nears  a solution  and  h'  (cu)  becomes  unstable. 

A tail  routine  such  as  Newton's  method  can  be  used  when  the  current 
point  is  near  an  equilibrium,  but  the  decision  as  to  when  to  institute 
this  tail  routine  is  by  no  means  an  easy  one.  Both  approaches  were 
experimented  with  and  the  homotopy  approach  was  found  to  be  much  more 
stable. 
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I 

IV . 3 . The  Path-FoIIo-wing  Algorithm 

Next  we  present  the  path  following  algorithm  i^ich  was  implemented 

i in  the  HRA  code  to  be  described  in  Chapter  V.  The  algorithm  is 

1 

f essentially  identical  to  that  presented  in  Section  we  include 

this  description  for  completeness.  Motivation  and  a verbal  explanation 
of  the  algorithm  are  given  in  Section 

Again  we  are  given  the  data  contained  in  a,  n,  v,  Gl,  G2,  c,  b, 

= (X  .t  ),  0.,  f as  defined  in  (l.9)^  and  (pd,  s)  describing  the 

0 (j.  V 0 

variable  currently  at  zero  which  is  allowed  to  increase.  Define 


bg(z,  0) 


t + b 

V 


G2 

I a(s) 


A + 
u 


c 

s 


if  pd  = 1 

if  pd  = -1  . 


Algorithm  3.1.  The  homotopy  retraction  algorithm  (HRA). 

0.  i = 0,  User  supplies  a and  e, . 

’ max  1 2 3 

1.  Calculate  F'(z^,0^)  = [H|h]  € . 

2.  Solve  Hy  = -h  (if  det  H = 0,  rearrange  columns  of 
[H|h]  -»  [H'  |h']  and  repeat  Step  2). 

Let  03  ;=  sgn(det  H) 
if  i > 0 go  to  Step  U. 

3.  8 :=  sgn((Vb^(z^,0^),  (y,l)>), 

03  6 *03,  go  so  Step  5. 


1 


t 
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k. 


5 :=  OD'O) 


5-  ''i ' ^ 

6.  The  tangential  approximation  to  F ^(O) 

T(a)  = a^a  + a^,  a e iR^  , 

where  a^  : = v 

and  a^  :=  (z^, e^). 

7.  Solve  the  n linear  equations  in  a defined  by 


(G1  , , T(a))  + b.  =0,  i e k 


1,  • 


(6.1) 


(G2  . , T(a))  + c.  = 0 , i e 1 

1, . i 


a,  , Q!  + a„  , ^ = 1 
l,m+l  2,m+l 

Let  a*  be  the  minimum  positive  real  root  of  all  these  equations. 
If  a primal  variaoie  hits  zero,  pd  : = 1 
If  a dual  variable  hits  zero,  pd  :=  -1 

Let  r ;=  the  index  in  (6.1)  corresponding  to  the  row  for  which 
T(cd^)  is  binding. 

Let  A :=  2. 

If  a,  a*  + a.  ^ = 1,  let  A : = 3. 

l,m+l  2,m+l 

Let  (z°,0°)  = T(a*) 

8.  If  or**-  < Q!  go  to  Step  9.  Otherwise,  let  o*  :=  a*/2,  (z*^,e^)  = 

T(o:*<-),  b^(z,fl)  5 (a^,  (z°,0^)  - (z,0)),  A ;=  1,  and  go  to  Step  10. 
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9.  If  pd  and  r describe  the  primal  variable  "p",  let  A ;=  3- 
If  A = 3 go  to  Step  13. 


10.  Define 


g(z,e) 


F(z,0) 


and  iterate 


= (z-^,0^)  - g'(z'^,0^)"^  S(z\e^),  I = 0,1,2,..., 


until  ||g(z^  )||  < Cg. 


11.  If  A = 2 go  to  Step  12. 

Otherwise,  let  i :=  i+1, 

(z^, 0^)  :=  (z'^'''^,0'^'''^),  and  go  to  Step  1. 

12.  Let  (z*,0*)  = (z^'''^,0'^^^) 

Check  that  Gl*t*  + b > -e,, 

V - y 

and  02 ‘>1*  + b > -e,  , 

u - 3 

and  0*  S [0,1]. 

If  so,  return  the  values  (2*^,0*,  pd,  and  r)  to  the  main  program. 
Otherwise,  let  pd  and  r define  the  most  infeasible  constraint, 


solve 


for  M e IR  , let 


b_(u-(z*,0*)  + (l-^)a  ) = 0 


(z°,0°)  = n(z*,0*)  + (l-n)a  , 


and  go  to  Step  10. 


13.  (Tail  Routine).  Perform  the  Newton  iteration, 

- f'(z^)"^  f(z^)  , I = 0,1,2,... 

until  ||f(z'^^^)||  < e^.  Let  (z*,0*)  :=  and  go  to  Step  12. 

Again,  several  remarks  are  in  order; 

a)  The  sizes  of  e,,  i € 3 ^ and  CK  may  be  chosen  to  suit 

' i’  max 

a particular  problem.  Typical  choices  mi^t  be  = 10“^,  Sg  = 

-U 

e,  = 10  , and  a =5.  The  choice  of  a will  radically  £rffect 

3 max  max 

algorithm  efficiency.  If  it  is  too  small,  many  Newton  iterations  will 
be  used  to  cross  larger  cells  when  it  would  probably  work  to  take  big 
Jumps  across  in  the  facet  which  the  curve  seem  to  hit.  If  a „ is 
too  large  when  F~^(0)  takes  a sharp  curve.  Step  10  may  not  converge. 

b)  (step  9)  "p"  is  the  variable  \Aich  the  original  auxiliary 

linear  program  was  maximizing.  Recall  from  the  discussion  in  Section 
V.2  that  p = ~ ^=1  ^1^“)  "near"  that  of 

the  theoretically  desirable  path  described  in  Section  IV, 1,  then 
f^(a3)  >0,  1 € m.  Thus,  if  the  path  hits  the  "p  = 0"  facet  we  are 
close  to  an  equilibrium  and  the  tail  routine  can  be  implemented. 

c)  It  is  worth  discussing  how  F'(z,  0)  Is  computed  in  St^  10 
assuming  that  a subroutine  is  available  for  computing  f'(7v^,t^).  From 
the  definition  of  F(z, 0)  in  (5.1)  it  is  apparent  that 


F'  (^  ,t  ,0) 

(i  V 


r 


The  calculation  is  not  much  worse  than  calculating  f*  (X  ,t  ) 

u V 

if  the  values  of  f(X  ,t  ) and  G2  A + c are  saved  from  the  calcu- 

u V u-  u u 

lation  of  F(X  ,t  ,0)  as  suggested  by  (5.1). 

u V 

d)  If  we  replace  the  first  line  in  Step  8 by 

"If  < a go  to  9,  otherwise,  let  cc*  :=  a _ , " 

TTlflX  QlcUC 

Then  this  algorithm  is  essentially  the  path  following  algorithm 
of  Chapter  I,  which  was  proven  to  be  convergent  given  e , and  a 

& fllElX 

are  small  enough.  This  course  of  action  would  result  in  an  inordinate 
amount  of  work  for  large  cells — work  which  is  not  necessary  in  most 
cases. 

IV, 4,  The  Cell  Switching  Decisions  and  Operations 

The  decisions  necessary  in  the  "main  program"  of  this  algorithm 
for  equilibrium  calculation  are  very  simple.  Initially  t^  = 0 and 
V = ^ so  any  one  of  the  first  m primal  variables  may  be  specified 
as  binding — we  choose  pd  = 1,  s = 1.  With  the  correct  initial  choice 
of  v^,  i £ all  budget  surpluses  will  be  positive.  An  argument 
similar  to  that  in  the  proof  of  Lemma  1.6  will  show  that  all  the  vari- 
ables t^,  1 £ V increase  Initially,  and,  hence,  the  initial  tangent, 
a^,  points  into  the  interior  of  the  initial  cell.  Thus,  no  problems 
are  encountered  due  to  the  fact  that  we  are  starting  out  at  the  vertex 
of  the  Initial  cell. 


L 
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Flowchart  for  cell  switching  and  basis  updating. 

(The  same  abbreviations  will  be  used  as  in  Section  VI. 4.) 


I. 


i 

* 

r 


START 


FIGURE  3.1 
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The  reader  is  directed  to  the  comments  following  Figure  VI. 5,1 


for  a discussion  of  the  abbreviated  description  of  the  algorithm  in 
the  boxes  of  Figure  5.1. 


CHAPTER  V 


RESULTS  FROM  COMPUTATIONAL  EXPERIMENTS 

In  this  chapter  we  report  the  results  of  numerical  experiments 
with  a code  which  implements  the  Bilinear  Complementarity  Algorithm 
(BCA)  of  Chapter  III  and  Homotopy  Retraction  Algorithm  (HRA)  of 
Chapter  IV.  A large  part  of  the  work  involved  in  this  experimentation 
was  in  the  construction  of  some  test  problems.  We  shall  discuss  how 
these  problems  were  generated  and  what  sort  of  preprocessing  was  necessary 
in  order  to  solve  the  problems  with  these  algorithms.  The  two  codes  are 
named  BCA  and  HRA,  respectively.  The  only  differences  in  the  codes  are 
those  absolutely  necessary  to  implement  the  different  path  definitions 
of  the  two  methods. 

Both  BCA  and  HRA  use  LPMl,  an  all-in-core  FORTRAN  linear  pro- 
gramming code  written  at  Stanford  by  J.  A.  Tomlin.  LPMl  stores  the 
problem  matrix  by  columns  packed  in  a vector  of  non-zeroes,  a vector  of 
the  same  dimension  giving  for  each  non-zero  coefficient  its  row  index, 
and  a vector  giving  for  each  column  the  position  of  its  first  nonzero 
element  in  each  of  the  two  above  vectors.  The  eta  file  is  stored  using 
the  same  principle.  It  uses  an  L-U  factorization  for  inverting  the 
basis  followed  by  product  form  updates. 
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L,  Computer  Implementation  of  the  Algorithms 

The  BCA  and  HM  codes  are  written  in  FORTRAN-TV  and  eire  compatible 
with  both  the  WATFIV  and  FORTRAN-H  compilers.  The  results  to  be  reported 
[ here  were  done  with  the  FORTRAN-H  compiler  with  full  optimization 

i 

[ (opt  = 2).  The  testing  was  performed  on  an  IBM  370-168  computer  located 

at  the  Stanford  Linear  Accelerator  Center.  BCA  consists  of  3,837  lines 
of  code  while  HRA  is  3,632  lines  long.  LPMl  occupies  1,819  lines  of 
each  of  these  programs. 

The  form  of  the  data  input  will  not  be  described  in  detail  here, 

but  most  of  the  data  is  input  as  the  auxiliary  linear  program  in  MPS 

stajidard  format.  A small  amount  of  additional  data  must  be  added  to 

describe  the  vector  C (*,!),  for  I = 1, ...,m,  which  contain  the 

th 

coefficients  necessary  to  calculate  the  budget  surplus  of  the  I con- 
sumer. The  details  of  input  and  output  are  contained  in  "A  Programmer's 
Guide  for  BCA  and  HRA;  Two  programs  for  computing  economic  equilibria," 
to  be  published  as  a technical  report. 

V,2.  Experimental  Design 

The  primary  objective  of  these  numerical  experiments  is  to 
demonstrate  that  a path-following  philosophy  along  with  an  exploitation 
of  the  linear-programming  structure  can  be  implemented  to  solve  some 
problems  of  moderate  size  (up  to  6 consumers  and  56  goods).  It  was 

hoped  that  the  run  times  would  be  short  enough  for  these  problems  that  i 
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it  would  be  clear  that  these  algorithms  show  promise  for  considerably 
larger  problems.  By  larger  problems  we  mean  those  concerned  with 
about  10  traders  and  250  goods.  A model  for  which  such  a capability 
would  be  useful  is  a model  of  international  trade  where  the  traders  are 
countries  or  groups  of  countries. 

A secondary  objective  is  to  compare  the  performance  of  the  BCA 
and  HHA  codes  as  they  process  the  same  set  of  test  problems.  The  only 
conclusion  we  can  make  is  that  the  BCA  performs  better  on  certain 
problems  and  worse  on  others,  in  comparison  with  the  HM  algorithm. 

One  of  the  biggest  difficulties  in  performing  this  experiment 
was  in  finding  or  manufacturing  suitable  test  problems.  Some  small 
examples  were  given  in  Wilson  [1976]  which  have  been  solved  by  hand. 

Two  of  these  were  solved  by  the  BCA  and  HEA  codes  mainly  to  verify  that 
the  codes  actually  were  computing  equilibria.  Three  problems  were 
adapted  from  the  three  equilibrium  problems  in  Scarf  [1973].  The  last 
two  problems  appear  for  the  first  time  here.  The  run  times  will  follow 
the  problem  descriptions. 

Problem  Descriptions 

Problem  1.  This  problem  is  due  to  Andreu  Mas-Colell.  There 
are  three  traders  and  two  commodities.  Each  trader  has  one  unit  of 
each  good  and  the  utility  functions  are,  respectively, 

min(x,  2y) 
min(2x,  y) 
min(4x,  5y) 
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where  x eind  y are  the  quantities  of  the  two  goods  consumed. 

One  can  easily  verify  that  the  demand  function  derived  from 
these  utility  functions  will  have  a ray  as  its  range.  Thus  the  activities 
of  the  traders  can  be  reduced  to  one  each  and  the  utility  functions  are 
reduced  to  linear  functions.  The  problem  data  is,  now 

= (1)  , i e 3 


g3  _ / *^3^ 
® ^ .20^ 


f 


= ( ^ ) ^ 


i € 3 . 


V*  = (1, 1, U)  is  determined  from  the  formula  (V.2.6).  We  choose  v^ 
to  be  slightly  less  than  vf  as  suggested  in  Section  V.2: 

= .9,  v^  = .95,  v^  = 3.92. 


This  example  is  used  to  show  that  an  equilibrium  problem  of 
this  form  can  have  rational  data  and  an  irrational  solution.  The 
equilibrium  prices  are  some  scalar  multiple  of  the  vector 

TT  = (1  + v/5,  1)  . 

The  prices  computed  by  both  algorithms  are  correct  to  six  decimal 
places, 

TT  = (.732051,  .2679^*9)  . 
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Problem  2.  This  example  is  due  to  Alan  Whisraan,  It  is  reported 
in  Wilson  [1976].  It  is  concerned  with  four  consumers  and  three  goods. 
The  utility  functions  are  linear,  but  the  matrices  are  non-trivial. 
Again,  there  are  no  firms. 

Problem  data; 


w = 2 


1111 

0 10  3 
12  3 0 
0 0 U 3 
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Problem  3.  This  is  the  first  of  three  examples  taken  from  Scarf's 


The  Computation  of  Economic  Equilibria  [1975].  This  is  a pure  exchange 
model,  while  the  next  two  include  production.  This  example  involves  10 
commodities  and  5 traders.  Each  trader  has  10  activities,  the  consumption 
I of  each  of  the  10  commodities,  so  the  activity  matrix  = I £ 

for  i = 1,  ...  , 5.  Each  consumer  has  an  initial  stock  w^  £ E of 
commodities  prior  to  trade  and  a utility  function. 

10  1/b.  1-l/b. 

i u. (y)  = sign(b  - !)•  Z a ^ y ^ , i £ 5 , 

j j=l 

i 

for  final  consumption.  This  is  a special  case  of  the  class  of  constant 

elasticity  of  substitution  or  CES  utility  functions.  For  convenience 

we  will  refer  to  this  functional  form  as  a CES  utility  function. 

Scarf's  algorithm  uses  the  demand  functions  derived  by  maximizing 

u.(y),  subject  to  E3*^7r.y.  < S^*^7r.w.  . to  derive  a mapping  vdiose  fixed 
1 J J -L  3 IJ 

point  corresponds  to  an  equilibrijm.  To  implement  either  of  our 
algorithms  we  must  calculate  a piecewise  linear  approximation  to  each 
consumer's  utility  function. 

From  an  econometric  point  of  view,  it  is  probably  Just  as  easy 
and  accurate  to  assess  a utility  function  in  a piecewise  linear  form 
as  in  some  standard  nonlinear  form.  However,  if  one  is  trying  to 
approximate  a nonlinear  functional  on  E by  linear  pieces,  it  may 
take  a very  large  number  of  pieces.  For  this  reason,  we  are  not  trying 
to  duplicate  the  answers  which  Scarf  computed,  although  in  some  cases 
we  come  close.  As  a rule  of  thumb,  we  compute  twice  as  many  pieces 
of  linearity  as  the  dimension  of  the  domain  of  the  utility  function. 
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In  this  example,  we  compute  PO  pieces  of  linearity  for  each  of  the  five 


consumers.  The  details  of  the  calculations  involved  in  this  prepro- 
cessing are  described  in  the  Appendix.  The  prices  obtained  in  the 
equilibrium  solution  were  not  very  close  to  those  computed  by  Scarf, 
but  that  is  to  be  expected. 

The  problem  data  is  given  below. 

Initial  Stock  of  Commodities 

Consumer 


1 

.6 

.2 

.2 

20 

.1 

2 

9 

5 

5 

15 

2 

.2 

11 

12 

15 

l4 

15 

l6 

5 

5 

9 

5 

.4 

9 

8 

7 

6 

5 

4 

5 

7 

12 

4 

1 

5 

5 

5 

5 

5 

5 

8 

5 

17 

5 

8 

1 

22 

10 

.5 

.9 

5.1 

.1 

6.2 

11 

Consumer 
1 1 

2 1 

5 9.7 

4 1 


Utility  parameters 


a.  . 
ij 


1 3 .1  .1  1.2  2 
1 111  1 1 
.15  .26  .28 
2 5 4 5 6 7 


5 1 15  11  9 4 .9  B 


1 1 .7 

111 
1 1 .2 
8 9 10 

1 2 10 


utility  parameters  b 


Consumer 


Problem  U.  This  is  the  Walrasian  equilibrium  model  on  pp.  109-113 
of  Scarf  [1973].  The  problem  involves  5 consumers,  6 commodities  and  a 
production  sector.  One  can  think  of  the  production  sector  as  being  a 
single  firm.  The  consumers  share  ownership  of  the  firm  in  the  sense 
that  they  own  the  factors  of  production — capital,  unskilled  labor,  and 
skilled  labor.  The  six  commodities  may  be  described  as  follows: 

1.  Capital  available  at  the  end  of  the  period. 

2.  Capital  available  at  the  beginning  of  the  period. 

3.  Skilled  labor. 
h.  Unskilled  labor. 

5.  Nondurable  consumer  goods. 

6,  Durable  consumer  goods. 

The  activity  matrix  for  the  firm  is 


-5.5  -5 


1.6  1.6  1.6 


-It  -5 


The  consumers  have  initial  endowments  of  commodities  2,  3,  4 
and  6 given  in  the  following  matrix: 

Commodity 


Consumer 

2 

3 

4 

6 

1 

3 

5 

.1 

1 

2 

.1 

.1 

7 

2 

3 

2 

6 

.1 

1.5 

4 

1 

.1 

8 

1 

5 

6 

.1 

.5 

.2 

The  consumer’s  utility  functions  are  of  the  same  form  as  in  the 
previous  problem.  The  specific  values  of  a. . are  given  in  the 
following  matrix; 


Commodity 


Consumer 

1 

2 

3 

4 

5 

6 

1 

4 

0 

.2 

0 

2 

3.2 

2 

0.4 

0 

0 

.6 

4 

1 

3 

2 

0 

.5 

0 

2 

1.5 

4 

5 

0 

0 

.2 

5 

4.5 

5 

5 

0 

0 

.2 

4 

2 

The  activity  matrix  is  composed  of  zeroes  except 

= 1 if  a^^  y'  0.  Of  course  zero  columns  can  be  omitted, 
elasticities  of  substitution  are  given  by 


that 

The 


b^  = (1.2,  1.6,  0.8,  0.5,  0.6) 


95 


Problem  5.  This  problem,  taken  from  pp.  113-119  of  Scarf 


[1975]  deals  with  four  consiimers  emd  fourteen  commodities; 

1.  Basic  agricultural  goods. 

2.  Processed  food. 

3.  Textiles. 

k.  Housing  services  and  heating. 

5.  Entertainment. 

6.  Housing,  end  of  period. 

7.  Other  capital,  end  of  period. 

8.  Steel. 

9.  Coal. 

10 . Lumber 

11.  Housing,  beginning  of  period. 

12.  Other  capital,  beginning  of  period. 

13.  Labor 

lU.  Foreign  exchange. 
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In  this  case  the  productive  sectors  involves  import  and  export 
activities.  The  activity  analysis  matrix  is  given  in  Scarf. 

The  utility  fixnctions  are  of  the  form 


7 a 

u.(y)  = n y, , , 


t=i 


'ij 


where 


7 

z 

i=l 


a. . = 1, 
13 


J = 1,  ...  , it 


These  Cobb-Douglas  utility  functions  are  homogeneous  of  degree  one, 

i.e.,  for  any  'h>Q,  Xu.(y)  = u.(>\y).  This  means  that  our  piecewise 

3 3 

linear  approximations  are  tangent  to  the  graph  of  u.  along  a ray 

3 

rather  than  at  a single  point.  As  one  might  expect,  this  means  that 
our  piecewise  linear  approximation  will  be  better  for  this  form  of 
utility  function  than  for  the  constant  elasticity  of  substitution 
utilities. 

The  coefficients  a.  . are  given  by 


Consumer 

1 

2 

3 

k 


1 

.1 

.2 

.5 

.1 


2 

.2 

.2 

.2 

.2 


Commodity 

5 ^ 

.1  .1 

.1  .1 

.3  .1 

.1  .1 


567 

.1  .3  .1 

.1  .1  .2 

.1  0.0  0.0 

.1  .1  .3 
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To  make  some  sort  of  comparison  between  the  Cobb-Douglas  and 


CES  (l)  utility  functions,  we  also  solved  this  equilibrium  problem 
using  the  former  type  of  utility  function.  We  used  intensities  a^^ 
proportional  to  those  in  the  table  above  and  let  b^  = 1.2,  i = 

Below  is  a comparison  of  the  equilibrium  prices  derived  from  this 
model  using  1)  Scarf's  results,  2)  our  results  using  PL  approximation 
to  the  Cobb-Douglas  utility  function,  and  5)  our  results  from  a PL 
approximation  to  a CES  utility  function.  (We  use  abbreviations  for 
the  names  of  the  commodities.) 


Equilibrium  Prices 


SCARF 

COBB- 

DOUGLAS 

CES 

1 

AGRl 

.0621 

.0613 

.0495 

2 

FOOD 

.0583 

.0595 

.0500 

3 

TXTL 

.0984 

.0977 

.0830 

k 

HSVH 

.0714 

.0706 

.0479 

5 

ENTR 

.0658 

.0650 

.0494 

6 

HEND 

.0624 

.0614 

.0674 

7 

CEND 

.0689 

.0722 

.0680 

8 

STEL 

.0981 

.0999 

.0831 

9 

COAL 

.0902 

.0888 

.0721 

10 

LUMB 

.0795 

.0773 

.0654 

11 

HBEG 

.0562 

.0552 

.1858 

12 

CBEG 

.0620 

.0650 

.0613 

13 

LABO 

.0365 

.0300 

.0328 

lU 

FEXC 

.0928 

.0956 

.0840 
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We  note  here  that  the  HEA  path  passed  through  only  3 cells  to 
reach  the  equilibrium  prices  in  column  2 while  it  passed  through  37 
cells  to  reach  the  equilibrium  prices  in  column  3-  The  computational 
results  reported  later  for  this  example  are  concerned  with  the  approxi- 
mation to  the  Cobb-Douglas  utility  fiinction. 

The  next  two  examples  are  new. 


Problem  6.  This  example  has  6 households  and  4 goods.  Five 
of  the  consimiers  have  piecewise  linear  utility  functions,  the  sixth  has 
a linear  utility  function.  Some  of  the  initial  endowments  of  the  con- 
siimers  are  zero.  In  this  case  the  entire  aiixiliary  linear  program 
matrix  will  be  given  to  show  how  the  piecewise  linear  utility  functions 
are  handled. 

One  notes  that  this  matrix  (Figure  2.1)  could  easily  be  put  into  a block 
diagonal  form  with  coupling  constraints,  but  since  our  LP  code  does 
not  tsike  advantage  of  such  structure  we  leave  it  in  the  form  such  that 
the  v^,  i = 1,...,  m are  in  the  first  m rows.  This  allows  us  to  consider 
the  slack  variables  t^,  i £ m as  the  first  m primal  variables  which 
eases  the  indexing  problems  in  the  code. 

The  initial  endowments  are  given  below. 

Commodity 


Consumer 

1 

2 

3 

k 

5 

6 


1 2 

■ 2 2 

1 )4 

0 0 

1 1 

1 2 

2 3 


3 

1 2 

3 2 

3 

5 0.5 

0 2 

2 1 


i 

I 


99 


-4.06 


r-l  OJ  -4-  ir\  VO 

>>>>>> 


VO  I— I iH 

VO  C7\  ON  VO 


f— 1 

rH 

rH 

LA 

CM 

rA 

fH 

rA 

1 

CM 

1 

CM 

'51 

W 

•V., 

^A 

rH 

rH 

CM 

I 


I 


CM 


I 

K^'l 


rH 

rH 

tT' 

CM 

CM 

A 

1 

rH 

rH 

rH 

rH 

rA 

1 

CVI 

rA 

rH 

rH 

.CM 

1 

H 

rH 

1 

CM 

f — 1 

h- 


w 


CM 


I— I rH 


CM 


LTN  CM 

l' 

CM 


CM 


r-l  CM  fO 

I I 


Ut  rTN  I-H 


l/N 

r-I 


CM 


LA 

H • 
I I 


CM  OJ  CJ 


<-l  CM  ro 

I 


O 

Al 

X 

.o' 

VI 

■p 

w 


X 

h5 

B 


o. 


100 


FIGURE  2.1 


We  will  call  this  problem  6a.  To  create  problem  6b  we  Just 
changed  three  entries  in  the  matrix  of  Figure  2.1. 

Problem  7.  This  is  a dynamic  model  in  which  three  consumers 
attempt  to  maximize  their  discounted  utility  subject  to  a budget  con- 
straint while  the  producers  maximize  profit.  The  model  is  an  adaptation 

« 

of  the  Mananaguay  model  of  W.  P.  Drews  [1976].  In  its  original  form 
it  was  a dynamic  linear  programming  model  of  the  economy  of  a developing 
country.  The  country  was  mythical  so  the  data  was  not  empirically 
supported,  but  the  numbers  were  reasonable  on  the  basis  of  past 
experience.  This  model  was  over  five  two-year  periods  and  included 
sophisticated  techniques  to  deal  with  distortions  due  to  the  finite 
planning  horizon.  The  model  offered  a rich  choice  of  objective  functions, 
including  one  which  maximized  a weighted  sum  of  consumption  and  one 
which  minimized  the  dependence  upon  imports. 

We  substantially  alter  the  model  1)  to  reduce  the  size  of 
the  problem  and  2)  to  make  the  expansion  in  the  consumption  sector 
which  an  equilibrium  formulation  allows.  We  reduced  the  size  of 
the  model  by  dropping  the  number  of  periods  to  three,  eliminating  import 
and  export  activities  and  simplifying  the  end-correction  mechanism. 
Arbitrary  levels  of  capital  stocks  are  defined  which  must  be  left  to 
the  society  at  the  end  of  the  third  period.  This  is  an  admittedly  poor 
method  for  handling  finite  horizon  problems,  but  its  simplicity  suits 
our  purpose. 
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We  will  not  present  the  data  of  the  problem  here,  just  a few 
statistics.  There  are  8 perishable  goods  in  each  period  (including 
three  kinds  of  labor)  and  8 capital  goods  in  each  of  the  three  periods 
plus  a fourth  set  of  8 capital  goods  which  is  bequeathed  at  the  end 
of  the  planning  horizon.  Since  the  same  good  in  two  different  time 
periods  is  considered  to  be  two  different  goods,  the  model  deals  with 
56  goods.  In  each  period  there  are  8 production  activities,  9 capital 
goods  construction  activities,  and  18  consumption  activities,  6 for 
each  consumer.  Hence,  there  is  a total  of  105  activities.  The 
utilities  are  discounted  at  a rate  of  15^  between  the  two  year  periods 
for  an  annual  rate  of  about  7.5^. 


Preprocessing 

Before  we  present  the  numerical  results  it  would  be  useful  to 
discuss  how  the  initial  utility  levels  v^,  i € m were  determined. 
For  the  pure  exchange  problems  1,  2,  3,  and  6,  we  used  the  procedure 
suggested  in  Section  V.2,  i.e.  solve  the  linear  programs 


V*  = maxfy^z^ |A^z^  < a^,  B^z^  < w^,  z^  > 0], 


i £ m , 


and  let  v^  = v*,  i G m.  In  general,  the  solution  of  these  m linear 
programs  took  about  as  much  computer  time  as  the  solution  of  the 
resulting  auxiliary  linear  program.  The  initial  utility  levels  v^^ 
derived  in  this  manner  were  very  satisfactory  for  these  models. 

For  the  models  with  production,  many  times  the  consumers  are  so  dependent 
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i 

I 

upon  the  other  consumers  and  the  productive  sector  that  it  does  not 
make  sense  to  determine  vt  as  suggested  in  Section  V.3.  In  these 
cases  we  either  guessed  at  the  initial  utility  levels  or  followed  the 
procedure  described  below; 

1,  Guess  at  the  equilibrium  prices  tt  of  the  consumers  initially  held 
commodities  w^. 

2,  Calculate  their  wealth  = ttw^. 

5.  Calculate  v^-  = max  u.  (z^) 

subject  to  TTZ^  < W^, 
using  the  nonlinear  utility  function  u^. 

This  value  vt  will  usually  be  less  than  the  equilibrium 
utility  level  for  consigner  i because  the  piecewise  linear  utility  over- 
estimates the  actual  utility  function  (see  Appendix  A).  If  a poor 
choice  of  initial  utility  level  causes  the  algorithm  to  fail  for  any 
reason,  a better  choice  can  usually  be  made  by  examining  the  output 
from  the  failure. 

Several  times  the  HRA  algorithm  failed  to  converge  because  of 
poor  initial  utility  levels.  The  BCA  failed  only  if  the  choice  of  some 
v^  was  so  large  that  f^(co)  = (C  ^,A)  - is  negative  after 

the  solution  of  the  auxiliary  linear  program.  We  will  report  one 
example  in  vrtiich  the  HRA  code  failed  but  the  BCA  code  worked. 
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V.3.  Numerical  Results 


1 


On  the  succeeding  pages,  tables  will  give  the  numerical  results. 

First,  though,  we  will  describe  the  various  measures  used  to  evaluate 
how  much  work  is  being  done  to  find  an  equilibrium  solution. 

1.  L.P.  iterations  and  L.P.  time:  The  solution  of  the 
axixiliary  linear  program  is  definitely  part  of  the  HRA  and  BCA  so 
the  amount  of  work  done  to  solve  it  must  be  measured.  If  the  utility 
levels  Vj^  are  very  good  the  largest  percentage  of  work  will  be  done 
in  the  LP  portion  of  the  algorithm.  The  LP  time  is  measured  using  the 
subroutine  LEFTIA  supplied  by  the  numerical  analysis  package  at  SLAC. 

The  accuracy  of  this  timer  seems  to  be  about  +.02  sec. 

i 

2.  Number  of  cells  traversed  and  path-following  time;  j 

The  number  of  cells  traversed  is  important  because  it  counts  the  number 

of  times  the  basis  and  the  superbasic  columns  are  updated.  This  work 
really  outweighs  the  work  done  to  follow  the  path  through  the  cell. 

The  path-following  time  measures  the  time  spent  after  the  auxiliary  LP 
is  computed  and  before  the  output  of  the  final  solution. 

3.  Scalar  function  calls:  This  refers  to  the  number  of  times 
any  of  the  consumers  budget  surplus  fimctions  are  evaluated.  In  the 
BCA  the  work  involved  is  on  the  order  of  the  dot  product  of  two  vectors 
in  IR  ^ (see  (111^4.9)),  where  d is  the  number  of  budget  surpluses 
currently  "at  zero."  In  the  HRA  there  are  always  m superbasic  vari- 
ables so,  on  the  aversige,  a little  more  work  is  done. 
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4.  Jacobian  evaluations:  To  compensate  for  the  changing 

dimension  of  the  BCA  here  we  counted  the  number  of  partial 

2 I 

derivatives  calculated  and  divide  by  m . In  the  HRA  we  just  co’ont  the  j 

number  of  times  a Jacobian  is  calculated.  Since  the  formulas  (III.4.U) 

j 

make  it  relatively  easy  to  compute  partial  derivatives  exactly,  that  is  - 

what  is  done.  Each  time  a jacobieui  is  calculated  a linear  system  is 
solved  in  either  code,  so  this  figure  represents  more  work  than  the 
evaluation  of  a Jacobian. 

5.  CPU  time:  This  includes  input  time,  which  is  not  negligible, 
and  the  output  time  of  the  equilibrium  solution  which  is  negligible. 

Thus  the  LP  time  8Uid  the  path-following  time  will  not  add  up  to  the 

CPU  time. 


Problem  Statistics 


Consumers, 


Problem 

Goods 

LP  size 

Density 

Iterations 

LP  time 

1 

3,2 

6x4 

50^ 

1 

.005 

sec 

2 

^,3 

8 X 11 

4o.9^ 

5 

.02 

sec 

3 

5,10 

ll6  X 56 

13.9^ 

95 

1.75 

sec 

4 

5,6 

52  X 36 

15.1^ 

33 

.54 

sec 

5 

4,l4 

75  X 57 

l4.6^ 

82 

.85 

sec 

6(a) 

6,4 

23  X 20 

18.5^ 

12 

.07 

se  . 

(b) 

16 

.08 

sec 

7(a) 

3,56 

97  X 107 

11.8^ 

96 

1.63 

sec 

(b) 

90 

1.63 

sec 

TABLE 

3.1 

103 


Comparison  of  Algorithms 
BCA  (Bilinear  Complementarity  Algorithm) 
HRA  (Horaotopy  Retraction  Algorithm) 


Problem 

Number 
of  Cells 
Traversed 

Path 

Following 

Time 

Scalar 

Function 

Calls 

Jacobian 

Evaluations 

CPU 

Time 

1 BCA 

1. 

*-» 

37 

4 

HRA 

2 

■H 

29 

2 BCA 

8 

.15 

45 

2 

.21 

HRA 

1 

,10 

16 

2 

.19 

3 BCA 

hk 

1.61 

l64 

7 

4.49 

HRA 

20 

.92 

319 

55 

3.90 

h BCA 

21 

.31 

82 

3 

.91 

HRA 

19 

.51 

387 

66 

1.17 

5 BCA 

i+6 

1,23 

166 

l4 

2.69 

HRA 

3 

.09 

44 

8 

1.44 

6a  BCA 

1)1 

.18 

194 

11 

.45 

HRA 

3 

.08 

90 

12 

.35 

6b  BCA 

Ik 

.21 

247 

l4 

.45 

HRA 

5 

.11 

130 

18 

CO 

7a  BCA 

4l 

1.83 

126 

13 

4.54 

HRA 

69 

3.47 

629 

173 

6.19 

7b  BCA 

33 

1.45 

70 

4 

4.18 

HRA 

60 

FAILED  TO 

COMPUTE  A SOLUTION^ 

TABLE  3.2 


^7b,  uses  the  same  data  as  7a.  A different  starting  point  caused  the 
HRA  to  diverge. 
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Comparison  of  HRA  on  Problem  5 with  Two  Different  Utility  Functions 


Number 
of  Cells 
Traversed 

Path 

Following 

Time 

Scalar 

Function 

Calls 

Jacobian 

Evaluations 

CPU 

Time 

Cobb- 

Douglas 

3 

.09  sec 

8 

m 

CES 

37 

l.Ul  sec 

559 

ll8 

TABLE  5,3 


] 

( 

I 

1 

The  solutions  calculated  by  the  codes  are  equilibrium  points  if 
they  satisfy  the  systems  of  equations  (ll,2.3)  and  (ll,2.3)  and  the  non- 
negativity of  the  variables.  Of  course,  no  finite  algorithm  can  calculate 
exact  equilibria,  in  general.  During  the  course  of  the  algorithm  the 

linear  inequalities  are  maintained  to  a tolerance  of  10  , i.e.,  in 

-4 

the  linear  approximation  Algorithm  III. 5.2,  = 10  . This  is  the 

default  tolerance  setting  in  LPMl  for  the  non-negativity  of  the  relative 
cost  coefficients.  The  final  application  of  Newton’s  method  for  calcu- 
lating a zero  of  the  budget  surplus  function  f(m)  terminates  when 
||f(a))l|  < 10  = 10  in  Algorithm  III. 5. 2).  These  tolerances 

are  among  the  user-supplied  parameters.  The  choice  of  is  rather 

arbitrary,  since  Newton's  routine  is  quadratically  convergent,  a smaller 

107 


4 


tolerance  will  not  appreciably  affect  run  time.  However,  since  the 


accuracy  of  the  veiriables  and  function  parameters  are  only  on  the  order 
-l4 

of  10  , it  would  not  be  appropriate  to  set  lower  than  that.  To 

set  Cg  smaller  than  10  may  appreciably  increase  the  run  time  of 
the  simplex  method  for  solving  the  auxiliary  linear  program.  Another 
tolerance  is  of  importance  in  these  algorithms.  The  termination 
criterion  for  each  Newton  routine  for  returning  to  the  curve  is  = lO’ 
Making  this  parameter  larger  may  reduce  the  run  time  marginally,  but 
one  runs  a risk  of  moving  outside  of  the  radius  of  convergence  surround- 
ing the  curve  with  the  next  step  along  the  tangent. 


V.  Conclusions 

The  conclusions  that  can  be  drawn  from  such  a limited  testing  must 
be  tentative  at  best.  However,  it  has  been  shown  that  some  medium-sized 
problems  can  be  solved  quite  quickly  using  path  methods.  In  our  opinion 
this  experience  shows  that  these  methods  show  promise  for  the  solution 
of  problems  with  up  to  10  consumers  and  200  goods.  This  is  only  specula- 
tion, though,  at  this  point. 

Referring  to  Table  3.2,  we  see  that  the  HRA  was  faster  6 times 
and  the  BCA  was  faster  on  3 problems.  This  indicates  that  there  can 
be  no  conclusion  drawn  concerning  which  method  is  better.  The  HRA 
appears  to  be  more  erratic.  Con^are  the  results  from  problem  5 and  7 
for  example.  The  increase  in  CPU  time  for  the  BCA  is  more  in  proportion 
to  the  change  in  problem  size  than  the  increase  for  the  HRA. 
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It  may  be  that  the  KRA  is  superior  for  problems  with  Cobb- 
Douglas  utility  functions.  We  cite  the  evidence  from  Problem  5 in 
Table  5.2  and  Table  3.3.  However,  more  testing  must  be  done  with 


problems  having  Cobb-Douglas  utility  functions  before  this  statement 
can  be  made  with  any  confidence. 

An  interesting  fact  not  reported  in  the  tables  above  is  that 
for  the  BCA,  the  number  of  bilinear  equations  vrtiich  are  binding 
increases  throughout  the  course  of  the  computation  on  such  problem. 

The  theory  does  not  guarantee  that  this  type  of  monotone  improvement 
occurs,  but  it  is  fortunate  that  this  behavior  seems  to  be  comnon. 

One  way  of  measuring  the  progress  of  the  HRA  is  to  note  how  B increases 
from  0^  e (0,l)  to  = 1 at  an  equilibrium.  In  several  examples 
6 increased  and  then  decreased  for  several  cells.  In  problem  7b, 

0^  = .256  increased  to  .697  before  decreasing  to  0,  causing  the  program 
to  stop  execution. 

The  HRA  passed  through  fewer  cells  than  the  BCA  in  every  problem 
but  Number  7j  this  advantage  must  be  large  in  order  for  the  HRA  to  take 
less  path  following  time  because  the  BCA  is  concerned  with  cells  of 
smaller  dimension  than  ra  most  of  the  time.  This  fact  explains  the 
relatively  small  number  of  scalar  function  calls  and  Jacobian 
evaluations  for  the  BCA  in  comparison  with  the  number  of  cells  traversed. 

This  comparison  of  algorithms  is  not  extensive  enough  to  allow 
us  to  make  authoritative  conclusions.  It  does  allow  us  to  say  that 
these  algorithms  show  promise  for  the  goal  of  solving  economically 
significant  equilibritun  models. 
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CHAPTER  VI 


SUGGESTIONS  FOR  FUTURE  RESEARCH 

We  close  this  work  with  some  ideas  for  continuing  this  investi- 
gation. Some  of  these  hypotheses  are  more  speculative  than  others 
but  they  all  deserve  to  be  mentioned. 

1.  Chapter  IV  gives  a convergence  proof  for  the 
Homotopy  Retraction  Algorithm  (HRA)  when  the  starting  point  >P  is 
the  zero  vector.  For  computational  efficiency,  we  solve  an  auxiliary 
linear  program  which  results  in  a starting  point  which  is  closer  to 
equilibrium  (Section  VJ , but,  as  problem  7h  indicates,  convergence 
is  no  longer  guaranteed.  Another  difficulty  is  in  choosing  the 
initial  utility  levels  v^,  i € m in  such  a way  that  a)  all  dual 
multipliers  i € m,  are  positive,  and  b)  all  budget  surpluses 

fj^(a)),  i e m are  positive. 

Although  the  author  has  tried  unsuccessfully  to  merge  other 
homotopy  deformations  with  the  piecewise  linear  economy  to  define  a 
convergent  path  method,  there  still  may  be  promise  in  such  an  enter- 
prise. It  is  certainly  possible  that  some  other  deformation  could 
define  a path  which  would  have  better  computational  properties  them 
the  HRA.  One  simple  alternative  would  be  to  involve  the  t variables 
(slacks  corresponding  to  the  utility  levels)  rather  than  the  A-variables 
in  the  deformation  F(a5, 0),  i.e., 
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F(a),  0)  = yf(co)  - (l-^)(t(a))  - t^)  . 

Other  possible  path  methods  which  could  be  adapted  to  this  problem 
are  the  strong  path  method  (Section  IV. 5,  Part  1)  or  any  of  the  class 
of  path  methods  discussed  in  Section  IV. 7 (Part  1). 

One  possible  remedy  for  the  initial  point  problems  of  the 
HRA  would  be  to  solve  the  auxiliary  linear  program  (11.2,3^  with  rather 
conservative  initial  utility  levels,  and  then  parametrically  increase 
those  v^  for  which  = 0.  One  would  continue  this  increase  until 
A^  was  positive  and  stop  before  f^(a))  went  to  zero  (A^  = 0 =»f^  > 0). 
When  these  conditions  were  satisfied  for  i = 1,  2,  . . . , m,  the  HRA 
would  be  instituted. 

2.  We  discussed  earlier  why  it  was  reasonable  to  consider 
models  with  few  consumers  and  piecewise  linear  utility  functions, 
but  interesting  extensions  to  our  theory  are  possible  which  would 
deal  with  the  problems  of  many  consumers  and  nonlinear  utility 
functions. 

In  competitive  economies  the  various  agents  make  decisions 
independently,  using  only  a little  common  information  (prices). 

Perhaps  some  sort  of  decomposition  principle  could  be  used  to  account 
for  the  effects  that  the  decisions  made  by  a large  number  of  households 
liave  on  the  system.  Another  method  to  ease  the  computational  burden 
would  be  to  take  advantage  of  the  block-diagonal  structure  of  the 
linear  program  by  using  a GGUB  algorithm  (see  Winkler  [1972]  for  a 
survey  and  unification  of  these  methods). 


; 

1 
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If  one  is  tiying  to  approximate  a nonlinear  utility  function 


in  many  variables,  better  methods  may  be  available  than  using  a 
globally  defined  piecewise  linear  function.  If  one  could  deal  with 
only  the  hyperplanes  which  are  binding  or  near  binding,  the  dimension 
of  the  problem  could  be  reduced  considerably.  The  rows  corresponding 
to  these  hyperplanes  are  columns  in  the  dual  auxiliary  program.  Perhaps 
some  sort  of  column  generation  algorithm  could  be  used  in  the  dual 
system  to  bring  in  a column  whenever  the  current  point  is  in  an  inaccurate 
portion  of  the  current  PL  approximation. 

3-  Nothing  has  been  said  in  this  dissertation  with  reference 
to  the  problem  of  stability  in  equilibria  theory  and  computation. 

If  the  economy  is,  for  some  reason,  displaced  from  equilibrium,  viiat 
process  will  bring  it  back  to  the  desired  point?  It  appears  as  though 
a path  of  the  form  of  the  horaotopy  retraction  algorithm  could  be 
defined  which  would  lead  to  the  equilibrium  for  any  point  in  a neighbor- 
hood of  the  equilibrium.  A sensitivity  analysis  of  the  auxiliary 
linear  program  could  yield  insight  into  the  stability  of  the  equilibrium. 


IIP 
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APPENDIX 

THE  CALCULATION  OF  PIECEWISE  LINEAR  UTILITY  FUNCTIONS 

There  are  clearly  many  possible  ways  to  compute  a piecewise 
linear  approximation  to  a given  utility  function.  We  will  present 
one  method  which  is  relatively  simple  and  easy  to  refine.  We  will 
not  give  any  detailed  justification  for  this  particular  method. 

The  epograph  of  a concave  utility  function  u: I? “ IR  is 
defined  as 

epo(u)  = f (x,t)  |x  € ]R  t e ]R  \ t < u(x)}  . 

One  could  evaluate  u(x^)  by  finding 

i 

t = sup  t = u(x  ) . 

(x*^,t)  € epo(u) 

We  will  define  a polyhedral  approximation  epo(u)  to  epo(u)  and  the 
piecewise  linear  approximation  u of  u can  be  evaluated  as 

fi(x^)  = sup  t . 

(x°,t) € epo(u) 

We  will  find  an  outer  approximation  using  hyperplanes  tangent  to  the 
boundaury  of  epo  u.  Assume  that  u is  differentiable.  Suppose 
t*^  = u(x®).  Then  the  tangent  hyperplane  to  epo(u)  = {(x,t)|t  - u(x)  < 0) 
at  (x®,t^)  is 

{(x,t)|<(-^(x°),l),  (x,t)  - (x°,t°))  =0)  , 
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or  letting  = t*^  - (^yli(x^)  ,x^)  we  have  that 

epo(u)  c:  {(x,t)|t  - ('7u(x°),x)  < c°}  . 

epo(u)  can  be  calculated  by  choosing  a number  of  points  of  tangeny 
{x^,  i =0,1,...,  m)  and  calculating  the  appropriate  gradients: 

epo(u)  = [(x,t)|t  - <Vu(x^),x)  < c^,  i - 0,...,m]  . (l) 

The  problem  which  remains  is  to  meike  some  sort  of  choice  of 
points  of  tangency.  Pick  two  utility  levels  and  which 

are  likely  to  be  in  the  range  of  those  under  consideration  in  the 
problem.  Given  a vector  v € = [x  € = 1,  x > 0),  it  is 

usually  quite  easy  to  find  scalars  such  that  u(k^v)  = U^.  Our 

procedure,  then,  is  to  choose  n vectors  v^,  . . . , v’^  scattered  in 
some  manner  in  the  unit  simplex  and  to  find  the  corresponding  points 
x^  = k^v^  such  that  u(x^)  = U^,  for  i £ n.  Then  we  choose  n 
more  vectors  ...  , v^'^  which  yield  x^  such  that  u(x^)  = U^, 

i £ 2n\  n.  The  points  x^,  i £ ^ will  be  the  points  of  tangency 
in  the  definition  of  epo(u),  (l). 

One  way  to  choose  v^,  i £ 2n  is  to  distribute  them  around  a 
point  V*  which  maximizes  the  utility  subject  to  a budget  constraint 
assuming  all  prices  are  equal,  i.e., 

V*  solves  maximize  u(v) 

n 

subject  to  E V,  < 1. 


X n 

We  arbitrarily  choose  V^,  . . . , v”  to  be  convex  combinations 
of  v^  and  each  of  the  vertices  of  If  e^  is  the  i^^  unit 

vector  for  some  i £ n,  then 

v^  = Qv*  + (l-a)e^  , i £ n . 

Let  d = (l/(n-l))e  - (l/(n-l))e^  where  e £ E is  a vector  of  ones. 
Then  v”^^,  is  chosen  to  be  some  convex  combination  of  v^  and  d^, 
the  centroid  of  S?  = fx  £ S*^|x.  = 0).  Let 

v”"^^  = QV*  + (l-a)d^  , i £ n. 


We  chose  a.  to  be  l/2,  but  other  choices  are  certainly  possible. 

To  illustrate  what  our  piecemse  linear  approximations  look  like 
we  must  use  examples  in  three  different  dimensions. 

The  first  illustrates  why  we  choose  points  of  tangency  at  two 
different  utility  levels  (Fig.  A.l);  in  general,  concave  utilities 
have  decreasing  marginal  returns. 


FIGURE  A.l 


FIGURE  A. 3 

It  may  be  that  the  hyperplanes  generated  with  this  particular 
choice  of  U^,  Ug,  and  a yield  answers  which  seem  to  be  bad  in  some 
sense.  In  that  case,  information  gained  from  the  solution  of  that 
problem  can  be  used  to  adjust  the  parameters  v*,  U^,  and  a so 

as  to  achieve  a better  approximation  to  u in  the  region  of  interest. 
With  this  refinement  the  equilibrium  problem  can  be  solved  again. 


117 


BIBLIOGRAPHY 


Arrow,  K.J.  and  G.  Debreu  [1954],  "Existence  of  an  Equilibrium  for  a 
Competitive  Economy,"  Econometrica, 22.  pp.  ?65-90. 

Arrow,  K.  and  F.  Hahn  [1971].  General  Competitive  Analysis.  San  Francisco: 
Holden  Day  Publishing  Company. 

Ari’ow,  K.J.  and  L.  Hurwicz  [195R],  "On  the  Stability  of  the  Competitive 
Equilibrium,  I",  Econometrica , 26,  pp.  522-52. 

Avila,  J.H.  [1974],  "The  Feasibility  of  Continuation  Methods  for  Nonlinear 
Equations,'  SIAM  J.  Numer.  Anal..  11,  pp.  102-22. 

Dantzig,  G.B.  [I963],  Linear  Programming  and  Extensions.  Princeton, 

New  Jersey:  Princeton  University  Press. 

Dantzig,  G.B. , _B.C.  Eaves,  and  D.  Gale  [1976],  "An  Algorithm  for  a 
Piecewise  Linear  Model  of  Trade  and  Production  with  Negative 
Prices  and  Bankruptcy,"  Department  of  Operations  Research, 

Stanford  University,  August,  I976. 

Debreu,  G.  [1951],  "The  Coefficient  of  Resource  Utilization," 

Econometrica.  19,  pp.  273-292. 

Debreu,  G.  [1959],  Theory  of  Value,  New  York:  Wiley. 

Drews,  W.P.  [1976],  "Demonstration  Examples  of  Trade-Payments  Model," 
iinpublished  communication. 

Eaves,  B.  C.  [1971b],  "The  Linear  Complementarity  Problem,"  Management 
Sciences,  17,  pp.  612-34. 

Eaves;,  B.C.  [1975],  "A  Finite  Algorithm  for  the  Linear  Exchange  Model," 
Technical  Report  SOL  75-4,  Stanford  University. 

Elken,  T.  [1977],  "On  the  Solution  of  Nonlinear  Equations  by  Path 

Methods,"  SOL  Technical  Report  77-25,  Department  of  Operations 
Research,  Stanford  University. 

Ginsburgh,  V.  and  J.  Waelbroeck  [1975],  "A  General  Equilibrium  Model  of 
World  Trade,  Part  I,  Full  Format  Computation  of  Economic 
Equilibria,"  Cowles  Foundation  Discussion  Paper  No.  4l2, 

Yale  University. 

Ginsburgh,  V.  and  J.  Waelbroeck  [1975],  "A  General  Equilibrium  Model 
of  World  Trade,  Part  II,  the  Empirical  Specification,"  Cowles 
Foundation  Discussion  Paper  No,  4l3,  Yale  University. 

Lemke,  C.E.  [1965],  "Bimatrix  Equilibrium  Points  and  Mathematical 
Programming,"  Management  Sciences,  11,  pp.  68I-89. 

Mantel,  R.  [1968],  "Toward  a Constructive  Proof  of  the  Existence  of 
Equilibrium  in  a Competitive  Economy,"  Yale  Economic  Essays, 

8,  pp.  155-196. 


118 


Murtagh,  B.A.  and  M.A.  Saunders  [1977%  "A  Large-Scale  Nonlinear 

Programming  System,  User's  Guide,"  Technical  Report  SOL  77-9, 
Department  of  Operations  Research,  Steinford  University. 

Ortega,  J.M,  and  W.C.  Rheinboldt  [1970],  Iterative  Solution  of  Nonlinear 
Equations  in  Several  Variables,  New  York:  Academic  Press. 

Quirk,  J.R.  and  R.  Saposnik  [I96S],  Introduction  to  Equilibrium  Theory 
and  Welfare  Economics,  New  York:  McGraw-Hill. 

Scarf,  H.  [I96O],  "Some  Examples  of  the  Global  Instability  of  the 
Competitive  Equilibrium,"  International  Economic  Review,  1, 
pp.  157-72. 

Scarf,  H.  [1967a],  "The  Approximation  of  Fixed  Points  of  a Continuous 
Mapping,"  SIAM  J.  Applied  Math. , I5,  pp.  I528-U3. 

Scarf,  H.  [1967b],  "On  the  Computation  of  Equilibrium  Prices,"  in  Ten 
Essays  in  Honor  of  Irving  Fischer,  ed.  Fellner,  et  al. , 
pp.  207-250,  New  York:  John  Wiley  and  Sons. 

Scarf,  H.  [1975]^  The  Computation  of  Economic  Equilibria,  New  Haven: 
Yale  University  Press. 

Smale,  S.  [I976],  "A  Convergent  Process  of  Price  Adjustment  and  Global 
Newton  Methods,"  Journal  of  Math.  Econ. , 3,  pp.  107-120. 

Wald,  A.  [19357>  "Uber  die  eindeutige  positive  Losbarkeit  der  neuen 
Produktionsgleichungen."  Ergibnisse  eines  Mathematischen 
Kolloquiums,  6,  pp.  12-20. 

Walras,  L.  [I87L],  Elements  d'Economie  Politique  Pure,  Lausanne:  Corbaz. 
[Translated  as  Elements  of  Pure  Economics,  trans.  W.  Jaffe, 
London:  Allen  and  Unwin,  195^. ] 

Wilson,  R.  [1976],  "The  Bilinear  Complementarity  Problem  and  Competitive 
Equilibria  of  Linear  Economic  Models,"  Technical  Report  SOL  76-2, 
Department  of  Operations  Research,  Stanford  University.  To 
appear  in  Sconometrica. 

Winkler,  C.  [197^],  "Basis  P’actorization  for  Block-Angular  Linear 

Programs:  Unified  Theory  of  Partitioning  and  Decomposition 
Using  the' Simplex  Method,"  SOL  Technical  Report  7^-I9> 

Department  of  Operations  Research,  Stanford  University. 


119 


UNCLASSIFIED 

SeCUSITV  CLASSIflCoTION  OF  THt*  P*Oe  n«»  tnltn*) 

REPORT  DOCUMENTATION  PAGE 

KP.AO  IN-STRl'CTIONS 

UKK  -KF,  COMHLETrNC,  KjK.M 

' RfiPORT  NLM.t"  [».  OOVT  ACCtSSlOH  NO. 

SOL  77-26  1 

i recip;Cn*"s  Catalog  NUMBtR 

« title  CAnd  SuS»/i/.> 

THE  COMPUTATION  OF  ECONOMIC  EQUILIBRIA  BY 

PATH  METHODS 

5 1 vrt  or  Rt-RORT  A RERIOO  COVEHED 

Technical  report 

• rcrformino  orc.  mcrort  number 

7.  AUTMORf.; 

Thomas  R.  Elken 

N00014-75-C-0267  • 
N00014-75-C-0865 

t PERAORM'NC  ON'.ASIZATICN  IIANI  AND  AODNESi 

Department  of  Operations  Research,  SOL 

1 Stanford  University 
j Stanford,  CA  94305 

«0.  RROORAU  ELEMENT.  RROJECT.  TASK 
AREA  A »ORK  UNIT  NUMBERS 

NR-047-064 

NR-047-143 

i't  CONTRoLLI.'.O  OFFICE  NAW..  ANO  ADDRESS 

j Operations  Research  Program 
! Office  of  Naval  Research 

1 Arlington,  VA  22217 

n.  REPORT  PATE  I 

October  1977 

tS.  NUMBER  or  RADES 

119 

t 14  monitoring  AGE'JCV  NAM*:  A ADOflCSS/lf  dllforonl  Inm  Controlling  OfHcm) 

1 

IS.  security  class.  (oI  thlt  fupotl) 

1S«.  OCClASSIPICATION/OOWNORAOinO 
SCHISOULC 

L- 

P*.  D'STPIBVJ  ttON  STAT^KMT  (<yt  thU  K^port) 


! This  document  has  been  approved  for  public  release  and  sale; 
} its  distribution  is  unlimited. 


• 17.  OISTRieuTtON  STATEMENT  fo/  »bstfsct  0fit0r0d  in  Block  70,  il  fram  Rmpoft) 


. I 

<_ - 4 

» IS.  supplementary  notes  I 

I i 


19.  KEY  WORDS  (Continxto  on  oi^m  H noeoooafy  m4  by  kioek  ntmbot) 

ECONOMIC  EQUILIBIRA  HOMOTOPY  RETRACTION  METHOD 

PIECEWISE  LINEAR  ECONOMY  PIECEWISE  LINEAR  APPROXIMATION  OF 

BILINEAR  COMPLEMENTARITY  PROBLEM  CONCAVE  UTILITY 
PATH-FOLLOWING  ALGORITHMS  COMPUTATIONAL  EXPERIMENTS 


\ 70.  abstract  (Continuo  on  rav*r»«  II  nmoooomy  on4  l^ontllf  by  bloob  mm^of) 


SEE  ATTACHED 


UNCLASSIFIED 

ItCUNiTv  CLAiSiriCATidM  or  thi*  paos  ONma  Bma 


SOL  77-26:  THE  COMPUTATION  OF  ECONOMIC  EQUILIBRIA  BY  PATH  METHODS, 
bv  Thom;is  R.  El  ken 


— Introduction  to  the  economic  equilibrium  model  is  given 
and  it  is  demonstrated  that  a path  method  can  be  used  to  compute 
equilibria  for  pure  exchange  economies  in  a nonlinear  setting. 

I Next,  a model  is  described  for  an  economy  in  which  the  utility 

I 

I functions  are  piecewise  linear  and  the  consumption  and  production  sets 
' are  polyhedral.  It  is  shown  that  an  equilibrium  for  this  economy  is 

I the  solution  to  a system  of  bilinear  equations  subject  to  certain 

5 

I linear  inequality  and  complementarity  constraints. 

1 Two  approaches  are  discussed  for  computing  equilibria  for  such 

I 

J economies.  The  first  is  the  bilinear  complementarity  algorithm  (BCA) 

' and  the  second  is  the  homotopy  retraction  algorithm  (HRA).  Convergence 
1 proofs  are  given  for  both  methods  using  the  general  theory  for  path 
j methods  described  above. 

The  BCA  and  HRA  have  been  implemented  as  computer  programs. 

Detailed  descriptions  of  the  algorithms  are  given,  and  the  results  of 

1 some  numerical  experiments  are  reported.  Seven  small  problems  were 
i 

j solved  by  both  algorithms No  conclusion  could  be  drawn  as  to  which 
' ^ algorithm  was  superior,  but  both  performed  well  enough  that  it  appears 
j that  much  larger  equilibrium  problems  also  can  be  solved  efficiently 
j by  these  methods.  
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